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Abstract 

A  new  method  for  approximating  anisotropic,  multi-group  scatter  cross  sections 
for  use  in  discretized  and  Monte  Carlo  multi-group  neutron  transport  is  presented.  The 
new  method  eliminates  unphysical  artifacts  such  as  negative  group  scatter  cross 
sections  and  falsely  positive  cross  sections.  Additionally,  when  combined  with  the 
discrete  elements  angular  quadrature  method,  the  new  cross  sections  eliminate  the  lack 
of  angular  support  in  the  discrete  ordinates  quadrature  method. 

The  new  method  generates  piecewise-average  group-to-group  scatter  cross 
sections.  The  accuracy  and  efficiency  for  calculating  the  discrete  elements  cross 
sections  has  improved  by  many  orders  of  magnitude  compared  to  DelGrande  and 
Mathews  (7)  previous  implementation.  The  new  cross  sections  have  extended  the 
discrete  elements  method  to  all  neutron-producing  representations  in  the  Evaluated 
Nuclear  Data  Files  (13). 

The  new  cross  section  method  has  been  validated  and  tested  with  the  cross 
section  generation  code,  NJOY  (13).  Results  of  transport  calculations  using  discrete 
elements,  discrete  ordinates,  and  Monte  Carlo  methods  for  two,  one-dimensional  slab 
geometry  problems  are  compared. 
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Efficient  and  Accurate  Computation  of  Non-negative 
Anisotropic  Group  Scattering  Cross  Sections  for 
Discrete  Ordinates  and  Monte  Carlo  Radiation  Transport 


I.  Introduction 

The  conventional  practice  for  evaluating  the  discretized,  multi-group  Boltzmann 
transport  equation  is  the  discrete  ordinates  angular  quadrature  method  with  truncated 
Legendre  expansions  representing  the  multi-group  cross  sections.  The  truncated 
Legendre  expansions  for  the  cross  sections  can  lead  to  negative  scatter  cross  sections. 
These  negative  scatter  cross  sections  can  give  rise  to  negative  scatter  sources,  which 
can  lead  to  inaccurate,  negative,  and  thereby  unphysical  angular  flux  solutions.  The 
unphysical  angular  flux  solutions  motivated  many  to  rely  on  comparatively  expensive 
Monte  Carlo  transport  simulations. 

Discrete  element  cross  sections,  recently  introduced  by  DelGrande  and  Mathews 
(7),  eliminate  the  negative  scatter  cross  section  artifacts,  but  efficient  and  accurate 
techniques  for  generating  these  discrete  elements  cross  sections  have  not  been  available. 

The  cross  sections  for  multi-group  Monte  Carlo  transport  also  inherit 
inaccuracies  by  attempting  to  reconstruct  non-negative  cross  sections  from  truncated 
Legendre  expansions  (5,  11). 

In  the  work  presented  here,  I  have  developed  algorithms,  implemented,  and  then 
validated  a  code  to  generate  multi-group  cross  sections.  The  code  improves  the 
efficiency  of  the  calculation  of  discrete  elements  cross  sections  for  discrete  ordinates 
transport  methods.  Furthermore,  it  provides  accurate  representations  suitable  for 
multi-group  Monte  Carlo  transport  methods.  The  utility  of  these  cross  sections  is 
demonstrated  using  various  discrete  ordinates  transport  calculations  with  two  test 
problems. 

1. 1 :  Multi-Group  Boltzmann  Transport  Equation 

The  multi-group  approximation  for  the  Boltzmann  transport  equation  (BTE)  is 


|VxS+s'("i’2(".W)=«'2(".W)+  a  0  gfe  (f  •  WxW)v  n  (r ,  W)  (l) 

g'E  1 

The  symbols  in  this  dissertation  are  summarized  in  the  table  of  symbols  in  the 
beginning  and  described  in  appendix  A.  The  notation  used  in  this  dissertation  follows 
the  conventions  of  Lewis  and  Miller  (7).  The  subscripts  in  the  multi-group  equation 
designate  the  group  ( g )  and  later  the  ordinate  or  element  ( n ).  The  spatial  dependence  is 
suppressed  in  all  further  equations.  The  superscript  designates  the  type  of  scatter  or 
the  total  cross  section.  For  further  development,  the  multi-group  scatter  cross  section 
is  defined  as 


s^(%xW)  = 


5  dEtfr  dE  ss  (E  ft®  F,%xW)f(£^) 

DE  ,  D Eq 

gt _ £ _ 

5  dEjX(E$) 


(2) 


where  F  (E  is  the  energy-dependent  spectral  weighting  function. 

The  multi-group  cross  section  can  be  further  discretized  in  angle  into  either  the 
discrete  ordinates  approximation  or  the  discrete  elements  approximation.  For  discrete 
elements,  discrete  ordinates,  or  multi-group  Monte  Carlo,  the  representation  of  the 
multi-group  cross  section  defined  in  equation  (2)  should  have  the  following  desirable 
attributes. 

1.  Non-negative,  but  not  strictly  positive 

2.  Account  for  all  incident  to  secondary  energy  group  pairs  for  all  scatter 
mechanisms 

3.  Efficiently  computed 

4.  Computed  within  user-specified  accuracy 

5.  Use  a  representation  that  is  flexible  enough  to  support  Monte  Carlo,  discrete 
ordinates,  and  discrete  elements  multi-group  transport  methods 


Multi-group  Monte  Carlo  transport  is  not  tested  in  this  work  and  is  presented 
clearly  in  the  reference  for  MCNP  (2).  The  discrete  elements  and  the  discrete  ordinates 
methods  are  introduced  and  compared. 
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1.2:  Discrete  Ordinates 


The  discrete  ordinates  method  samples  the  angular  flux  at  given  directions  and 
uses  a  quadrature  rule  to  integrate  the  angular  flux  to  get  the  scalar  flux. 

The  one-dimensional,  slab  geometry,  multi-group  discrete  ordinates  equation  is 


m  — y  +  sl  v„  „CO=  q„  „  +  a  a  s  a  ,,  V  . 

’  dx  1,8  n,g  U  U  nitugig"  n<kgi 

nl  gi 

where  the  ordinate-to-ordinate,  group-to-group  scatter  cross  section  is 

5  dEtfr  dEss(Et®  E,W^Wny(Efy 


DEst  DEs 
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Discrete  ordinates  conventionally  uses  the  Legendre  moments  of  the  cross 


sections,  s 


,  directly  by  computing  the  Legendre  moments  of  the  angular  flux,  f. , 


and  combining  those  with  the  cross  section  Legendre  moments  to  get  the  scattering 
source  for  each  ordinate,  n. 


The  series  of  approximate  discrete  ordinate,  angular  flux  solutions  to  the  multi¬ 
group  BTE  equation  do  not  converge  uniformly.  As  more  points  are  added  to  the 
angular  quadrature,  the  discontinuous  solution  (in  multi-dimensional  transport)  for  the 
angular  flux  suffers  from  a  Gibbs  phenomenon.  Therefore,  increasing  the  number  of 
ordinates  used  in  a  solution  does  not  guarantee  a  better  answer. 

The  discrete  ordinates  approximation  can  also  skip  energy  groups  in  the  down- 
scatter  arising  from  an  incomplete  angular  quadrature  approximation,  called  lack  of 
angular  support.  DelGrande  and  Mathews  (7)  gave  an  example  of  the  lack  of  angular 
support  for  discrete  ordinates  for  multi-group  scatter  cross  sections  that  skipped  109  of 
the  next  lower-energy  groups  out  of  175  groups.  Lack  of  angular  support  can  lead  to 
computational  artifacts  where  an  incorrect  scalar  flux  is  calculated. 


1.3:  Discrete  Elements 

The  discrete  elements  method  (7)  integrates  the  BTE  over  discrete  angular 
elements  with  a  piece-wise  constant  representation  of  the  angular  flux.  The  scalar  flux 
is  then  a  sum  over  each  of  the  angular  flux  elements.  The  resulting,  angularly 
discretized,  multi-group  BTE  has  the  form 
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DelGrande  and  Mathews  (7)  referred  to  the  cells  of  the  Cartesian  product  of  the 

energy  and  direction  meshes  as  bins.  Thus,  s  ,  ,  is  the  cross  section  for  scatter 

n  in, gig 

from  bin  ( n\g ')  to  bin  ( n,g ).  The  bin-to-bin  cross  section  is  defined  as 
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DelGrande  and  Mathews  (7)  showed  that  the  discrete  elements  approximation 
has  several  properties  that  are  advantageous  compared  to  the  discrete  ordinates 
approximation.  The  first  property  is  the  convergence  of  the  discrete  elements 
approximation  to  the  solution  of  the  spatially  discretized  problem  because  the 
approximation  is  equivalent  to  a  Riemann  integral  for  the  scalar  flux  in  the  limit  as 
N  ®  ¥  .  Therefore,  as  the  angular  mesh  is  refined,  the  refined  answer  is  more 
accurate  than  the  previous  answer.  Conversely,  the  discrete  ordinates  angular 
approximation,  being  based  on  the  spherical  harmonics  functions,  converges  only  point- 
wise  and  refinements  in  angle  for  discrete  ordinates  may  not  improve  the  scalar  flux 
solution. 

The  second  important  property  of  discrete  elements  is  complete  angular  support. 
Discrete  elements  have  non-zero  cross  sections  between  all  allowed  energy  groups  for 
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all  possible  scatters.  Again  conversely,  the  discrete  ordinates  approximation  need  not 
have  all  possible  energy  groups  included  in  the  cross  section  representation. 

1.4:  Motivation 

Multi-group,  non-negative,  anisotropic  scatter  cross  section  methods  for 
neutron  transport  have  previously  been  either  too  expensive  computationally  (7)  or  too 
restrictive  in  application  (2)  to  effectively  implement  for  real  materials  and  for  any 
energy  refinement. 

1.5:  Goal  of  the  Research 

The  goal  is  to  develop,  implement,  and  validate  efficient  and  accurate 
computational  methods  for  converting  the  Evaluated  Nuclear  Data  Files  version  B-VI 
(ENDF/B-VI)  (18)  to  multi-group  cross  sections,  without  an  intermediate  truncated 
Legendre  expansion,  in  forms  suitable  for  Monte  Carlo,  discrete  ordinates,  and  discrete 
elements  transport  calculations. 

1.6:  Scope 

I  include  the  reactions  in  the  ENDF/B-VI  data  where  the  neutrons  are  both 
incident  and  secondary  particles.  The  new  methods  for  calculating  multi-group 
anisotropic  scattering  cross  sections  are  validated  and  compared  to  other  methods.  The 
efficacy  of  these  cross  sections  for  discrete  ordinates  and  discrete  elements  one¬ 
dimensional  slab  geometry,  multi-group  transport  are  evaluated. 

1.7:  Assumptions  and  Limitations 

The  approximation  for  Doppler  broadening  of  the  cross  section  used  by  NJOY 
(13)  is  assumed  to  be  adequate  and  I  use  it  here.  Additionally,  the  R-Matrix,  hybrid  R- 
Matrix,  Adler-Adler,  and  Kalbach-Mann  representations  of  the  scatter  cross  sections 
that  are  used  for  some  isotopes  in  ENDF/B-VI  are  not  included.  They  are  not  required 
for  the  evaluation  and  demonstrations  presented  here,  and  are  left  for  a  future  effort. 
Utility  of  these  cross  sections  for  multi-group  Monte  Carlo  calculations  is  self-evident 
and  is  not  demonstrated  due  to  time  constraints.  Use  of  these  cross  sections  with  multi¬ 
dimensional  discrete  ordinates  and  discrete  elements  is  not  demonstrated,  also  due  to 
time  constraints. 
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1.8:  Approach 


This  work  starts  with  the  discrete  elements  approximation  introduced  by 
DelGrande  and  Mathews  (7).  They  used  Monte  Carlo  numerical  quadrature  to 
approximate  the  six  dimensional  integrals  for  the  bin-to-bin  cross  sections.  My 
approach  is  to  reduce  the  computational  cost  by  splitting  each  six  dimensional  integral 
in  equation  (9)  into  two,  three-dimensional  integrals  and  introducing  a  scattering  cross 
section  operator.  All  the  necessary  integrals  are  approximated  using  deterministic 
numerical  quadratures  to  improve  the  computational  efficiency  and  accuracy  as 
compared  to  the  Monte  Carlo  numerical  estimate.  User-set  tolerances  are  used  in  the 
deterministic  quadratures  to  ensure  the  desired  accuracy  is  achieved. 
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II.  Numerical  Approximations  for  Scatter  Cross  Sections 

This  chapter  presents  numerical  approximations  for  the  scatter  cross  section. 
The  conventional  truncated  Legendre  expansions  are  presented  first,  followed  by  the 
Monte  Carlo  technique  used  by  DelGrande  and  Mathews  (7)  to  estimate  bin-to-bin 
cross  sections.  A  scatter  cross  section  operator  is  introduced  and  is  shown  to  recover 
both  the  truncated  Legendre  expansion  and  the  bin-to-bin  cross  sections.  An  element- 
to-element  conditional  probability  is  introduced  for  the  discrete  elements  quadrature. 
Piecewise-average  group-to-group  scatter  cross  sections  are  introduced,  hereafter 
referred  to  as  PAX  cross  sections.  Finally,  the  failure  of  a  point-wise  method  for 
tabulating  the  group-to-group  cross  sections  is  presented. 

PAXk  is  a  symbol  denoting  piecewise-average  cross  sections  with  a  uniform 
mesh  of  K  (equal-width)  intervals  of  m,  ie  K  pieces.  In  a  complete,  formal  usage,  one 
could  refer  to  DE12/PAX61.  analogously  to  S12/P11.  But,  as  long  as  K  is  large  enough 
that  the  numerical  approximation  is  negligible,  DE12  should  suffice. 


II.  1 :  Truncated  Legendre  Expansions 


Current  practice  is  to  use  truncated  Legendre  expansions  of  the  group-to-group 
cross  sections  in  discretized  neutron  transport.  The  angular  domain  is  typically 
approximated  with  discrete  ordinates,  as  discussed  in  chapter  1.  These  two 
approximations  have  their  own  computational  artifacts  that  are  important  to 
distinguish.  The  coefficients  for  the  expansion,  or  the  Legendre  moments,  are 


5  dm  5  dEi  q  dE  ss  (Ei®  E,m) F  (E  ()P[  (ni) 


-  1  DE 


g* 


D  E 


gt 


1, gig 


5  dE&(Efy 

D  E  . 
gt 


(10) 


Having  tabulated  the  Legendre  moments  through  order  L,  the  group-to-group  cross 
section  can  be  recovered  (approximately)  as 


s  ,  (ni) » 


L 

o 

a 

1=0 


(21+  1), 


hg<g 


Pl  (m) 


(11) 


This  approach  has  the  advantage  of  compact  storage.  However,  because  the 
Legendre  polynomials  (for  />0)  are  not  non-negative,  the  recovered  cross  sections  can 
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be  negative  and  in  general  have  many  regions  of  negativity.  Simply  increasing  the 
order  of  the  expansion,  L,  does  not  necessarily  help,  and  may  even  exacerbate  the 
problem  in  certain  regions. 

To  show  some  problems  with  truncated  Legendre  expansions,  I  used  the  10B 
isotope  with  various  group-to-group  pairs  and  various  energy  group  structures.  The 
first  two  examples  use  the  30-energy  group  structure  given  in  appendix  D.2  to  show 
that,  in  practice,  it  is  not  possible  to  select  an  order,  L>0,  to  use  for  a  material  and  all 
energy  groups  pairs  and  maintain  positivity.  The  two  examples  have  either  negative 
regions  for  L  small,  but  are  strictly  positive  for  larger  L,  or  negative  regions  for  all 
orders  of  L>0.  All  of  the  expansions  for 

the  cross  sections,  from  P0  through  Pn,  are  shown  on  each  plot  and  compared  to  the 
PAX64  cross  section. 


Cosine  of  scatter  angle  (g) 

Figure  1:  Negative  regions  for  low  order,  L,  with  group  1  to  group  5 


Cosine  of  scatter  angle  (m) 

Figure  2:  Negative  regions  for  all  L>0  using  group  1  to  group  2 

These  examples  demonstrate  that  the  choice  of  the  truncation  order  is  either 
dependent  on  the  group-to-group  pair,  as  in  Figure  1,  or  cannot  be  chosen  to  guarantee 
non-negativity  for  L>0,  as  in  Figure  2. 

The  truncated  Legendre  expansion  through  order  L  will  often  develop  negative 
regions,  as  the  energy  group  structure  is  refined.  The  10B  isotope  for  group  1  to  group 
5,  shown  in  Figure  1,  is  an  example  of  the  onset  of  negative  regions  as  the  energy  group 
structure  is  refined.  The  example  in  Figure  1  favors  a  high-order  truncated  Legendre 
expansion  because  the  group-to-group  cross  section  is  strictly  positive  and,  indeed,  all 
expansions  greater  than  Pi  are  strictly  positive.  In  the  next  figures,  only  the  Pn 
Legendre  expansion  is  shown  and  compared  to  the  PAX64  cross  sections  for  a  refined 
energy  group  structure  using  117  energy  groups,  described  in  appendix  D.4.  The 
secondary  energy  group,  group  5,  has  been  refined  into  four  groups  of  equal  width  in 
lethargy,  with  the  incident  energy  group  held  constant. 
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Cosine  of  scatter  angle  (p) 

Figure  5:  10B  cross  section  for  group  1  to  16  in  1 17-group  structure 


-1  0  1 


Cosine  of  scatter  angle  (|l) 

Figure  6:  10B  cross  section  for  group  1  to  17  in  1 17-group  structure 

The  truncated  Legendre  expansions  may  perform  adequately  for  the  occasional 
combination  of  material,  Legendre  order,  and  energy  group  structure.  But  changing 


any  one  of  these  contributors  to  a  previously  adequate  combination,  for  example  the 
energy  group  in  the  figures  above,  can  lead  to  regions  of  negativity  in  the  truncate 
Legendre  expansion.  Therefore,  although  the  Legendre  moments  are  efficient  to 
generate  and  use  for  discretized  transport,  a  truncated  Legendre  expansion  does  not 
guarantee  the  most  important  quality  for  transport  data — non-negativity. 

Another  property  of  the  group-to-group  cross  sections  is  discontinuities  in  the 
first  derivative.  And  consequently,  the  globally  smooth  truncated  Legendre  expansions 
do  not  converge  uniformly.  This  point-wise  convergence  of  the  truncated  Legendre 
expansion  is  another  disadvantage  of  the  approximation  because  increasing  the  number 
of  moments  included  in  the  expansion  to  recover  the  group-to-group  cross  sections  does 
not  guarantee  a  more  accurate  solution  throughout  the  entire  domain  of  the  expansion. 

When  the  two  approximations,  discrete  ordinates  and  a  truncated  Legendre 
expansion  for  the  cross  section,  are  combined,  several  artifacts  are  obscured. 
Additionally,  the  combination  of  the  two  approximations,  neither  of  which  converges 
uniformly,  leads  to  a  certain  art  whereby  the  order  of  the  Legendre  truncation  and  the 
number  of  ordinates  to  use  is  divined  through  experience. 

The  lack  of  angular  support  in  the  discrete  ordinates  is  obscured  because  the 
truncated  Legendre  expansion  is  non-zero  for  all  points  between  -1  and  1  except  for  the 
finite  number  of  nodes.  The  combined  approximation  will  have  both  negative  cross 
sections  for  certain  directions  and  positive  cross  sections  for  other  directions  where  the 
group-to-group  cross  section  would  be  zero.  Thus,  the  lack  of  angular  support  is 
obscured. 

Combining  truncated  Legendre  expansions  and  discrete  ordinates  angular 
quadratures  can  lead  to  negative  scalar  fluxes  in  two  ways.  Using  the  equation  to 
calculate  the  scalar  flux 

N 

fg  =  a  wnyn,g >  (12) 

n  =  1 

having  either  negative  weights,  w«,  or  negative  angular  fluxes  can  lead  to  negative 
scalar  fluxes.  The  negative  angular  fluxes  can  arise  from  having  a  discrete  ordinates 
quadrature  set  with  an  angle  between  two  ordinates  that  falls  in  a  region  of  negativity 
in  the  truncated  Legendre  expansion  for  the  group-to-group  scatter  cross  section.  This 
ordinate-to-ordinate  pair  with  a  negative  value  for  the  scatter  cross  section  can  lead  to  a 
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negative  scatter  source.  If  no  other  source  is  present  in  the  spatial  cell,  then  the  angular 
flux  solution  will  be  negative.  Having  an  angular  quadrature  set  with  negative  weights 
(although  this  is  always  avoided  in  practice)  could  also  lead  to  negative  scalar  fluxes. 
Therefore,  I  reject  discrete  ordinates  with  truncated  Legendre  expansions  for  the  cross 
section  because  the  two  important  properties  of  the  BTE  can  be  lost  using  this 
combination  of  approximations — non-negative  fluxes  and  non-negative  input  data. 


II. 2:  Evaluation  of  Bin-to-Bin  Cross  Sections  by  Multidimensional  Monte 
Carlo  Numerical  Integration 


DelGrande  and  Mathews  (7)  used  Monte  Carlo  numerical  integration  to 
generate  the  bin-to-bin  cross  sections.  This  bin-to-bin  cross  section,  restated  from 
chapter  1,  is 


s 


n  in, gig 


Q  (/Wq  dW  q  dEtfr  dEss  (Ej®  £,%x\y)F(£$ 

DW  DW„  D  E  .  D  E„ 

n$ _ n_ _ gt_ _ £ _ 

DW  o  dEt?(E$) 


(13) 


This  integral  is  slow  to  converge  with  Monte  Carlo  and  the  stochastic 
convergence  does  not  guarantee  that  all  possible  bin-to-bin  scatters  are  taken  into 
account  without  using  an  infinite  number  of  samples.  And,  of  the  bin-to-bin  scatters 
that  are  taken  into  account,  many  bin-to-bin  pairs  will  have  large  stochastic  errors  due 
to  the  slow  convergence  of  Monte  Carlo  integration.  In  practice,  DelGrande  and 
Mathews  (7)  assumed  that  scatters  with  low  probability,  and  consequently  large 
stochastic  error,  would  not  contribute  significantly  to  the  transport  result.  Scatters 
with  low  probability  can  be  the  dominant  source  of  particles  at  some  energies  in  some 
locations  of  some  transport  problems.  I  present  an  example  of  this  in  a  shield 
penetration  problem  in  chapter  6.  I  reject  the  Monte  Carlo  evaluation  of  the  bin-to-bin 
cross  section  because  of  its  slow  convergence  and  large  stochastic  error. 

To  reduce  the  stochastic  error  and  improve  efficiency  and  accuracy,  variance 
reduction  techniques  were  investigated  for  use  in  the  Monte  Carlo  evaluation  of  the  bin- 
to-bin  cross  section.  As  the  various  ENDF/B-VI  scatter  mechanisms  were  investigated, 
several  of  the  integrations  had  either  analytic  solutions  or  could  be  deterministically 
integrated.  The  replacement  of  the  Monte  Carlo  integrations  using  deterministic 
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quadratures  or  analytic  integrations  continued  following  the  textbook  on  variance 

reduction  for  Monte  Carlo  by  Hammersley, 

It  should  almost  go  without  saying,  if  it  were  not  so  important  to  stress  it,  that 
whenever  in  the  Monte  Carlo  estimate  of  a  multiple  integral  we  are  able  to 
perform  part  of  the  integration  by  analytical  means,  that  part  should  be  so 
performed.  As  in  some  other  kinds  of  gambling,  it  pays  to  make  use  of  one’s 
knowledge  of  form.  (9) 


Eventually,  it  was  determined  that  all  of  the  nested  integrations  in  equation  (13)  could 
be  integrated  either  deterministically  or  analytically.  The  Monte  Carlo  quadrature  used 
by  DelGrande  and  Mathews  (7)  was  abandoned,  leading  to  the  methods  used  in  section 
II. 4. 

II. 3:  Scatter  Cross  Section  Operator 


The  two  methods,  bin-to-bin  cross  section  evaluation  with  Monte  Carlo 
integration  and  truncated  Legendre  expansions  for  cross  section,  can  be  combined  and 
extended  using  the  scatter  cross  section  operator  that  I  define  as 


q  dm  q  dE  q  dE  0>- s  (E  <t®  E,m) F  (E  f)g 

-  1  D  E  .  D  E„ 

_ gt _ £ _ 

5  dE$(E$) 


(14) 


where  g  is  any  real  function  of  m  defined  everywhere  in  the  interval  £-1,1]. 


S  ,  is  a 


mapping  from  the  function  space  which  comprises  its  domain  to  the  real  line.  It  is  a 
linear  functional  as  defined  by  Stakgold  (20). 


(m)  = 


5  i/W0q  dWd(m-  W0xW) 
DW  DW„ 

n  j- _ _ _ 

DW  , 

n 


(15) 


is  the  conditional  probability  that  a  particle,  uniformly  distributed  in  D  W  will  scatter 


into  DM(,  given  that  it  does  scatter  and  that  the  cosine  of  the  angle  of  that  scatter  is  m. 
Because  this  conditional  probability  is  entirely  determined  by  the  choice  of  a  partition  of 
the  unit  sphere  into  elements  of  solid  angle,  ,  ie.  by  the  discrete  elements  angular 
quadrature  set,  h  (in)  can  be  pre-computed  and  tabulated. 
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The  element-to-element  conditional  probability  function  has  been  numerically 
evaluated  and  tabulated  using  a  Monte  Carlo  numerical  integration.  It  has  also  been 
numerically  evaluated  using  a  Gauss-Chebyshev  numerical  quadrature  for  the  special 
case  of  one-dimensional  transport  (14). 

The  scatter  cross  section  operator  acting  on  the  function  h  ( m )  is 


5  dm  5  dE  (  q  dEss  (Et®  E,  m) F  (E$)  q  dM  q  dWd(m  -  %x\v) 
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Rearranging  the  integrals  and  performing  the  integration  over  the  delta  distribution 
function  introduced  in  equation  (15)  gives  the  equation 
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The  bin-to-bin  cross  section  is  therefore 

S  n  in, gig  ^ gig  in' 

The  Legendre  moments  for  the  truncated  Legendre  expansion  of  the  scatter 
cross  section  can  be  obtained  by  operating  on  P.  (in)  with  S: 
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which  is  the  equation  for  generating  the  Legendre  moments  of  the  scatter  cross  section. 

Piecewise-average  group-to-group  scatter  cross  sections  (PAX)  can  be  obtained 
using  the  scatter  cross  section  operator: 


,  =  S  , 

gig,k  g  ,g 


H(m-  mk_  ,  )//  (mk  -  m) 
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where  H  is  the  Heaviside  function,  k  =  OK  K ,  and 


mk  =  -  1  +  — .  (21 

PAX  cross  sections  have  several  advantages  compared  to  either  truncated 
Legendre  moments  or  Monte  Carlo  evaluation  of  bin-to-bin  cross  sections.  They  are 
non-negative  and  converge  uniformly.  The  PAX  cross  sections  can  be  used  for  discrete 
elements,  discrete  ordinates,  or  multi-group  Monte  Carlo  transport  having  been 
calculated  once  for  given:  material,  energy  group  structure,  F  (E  fy,  number  of  pieces, 

and  temperature.  In  the  remaining  chapters,  I  have  developed,  validated,  and 
demonstrated  an  algorithm  to  calculate  the  PAX  cross  sections  using  deterministic 
quadratures  for  the  integrations  in  equation  (20).  My  algorithms  use  deterministic 
quadratures  to  control  the  quadrature  error  introduced  in  approximating  equation  (20). 
PAX  cross  sections  are  accurate  for  each  group-to-group  pair  because  each  mechanism 
represented  in  ENDF/B-VI  and  each  group-to-group  pair  are  calculated  independently. 

II. 4:  Applications  of  Piecewise-Average  Group-to-Group  Scatter  Cross 
Sections 


II.4.1:  Bin-to-Bin  Cross  Sections 


PAX  cross  sections  can  be  used  in  a  very  efficient  way  to  approximate  the  bin- 
to-bin  cross  sections: 


where 
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(23) 


If  the  piecewise  conditional  scattering  probabilities 


n  j in  ,k 


are  pre-computed  and  stored 


for  a  particular  discrete  elements  angular  quadrature  set  and  choice  of  K,  and  similarly, 
the  PAX  cross  sections  T  ,  ,  are  pre-computed  and  stored,  then  the  set  ofj  , 

gtg,k  r  r  n  fin  ,g/ig 

values  is  obtained  by  the  simple  tensor  contraction  in  equation  (22)  with  no  redundant 
calculations. 
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II.4.2:  Legendre  Moments 


Likewise,  a  numerical  approximation  to  the  Legendre  moments  of  the  scatter 
cross  section  in  equation  (19)  using  the  PAX  cross  section  is 

o  (21  +  1) 
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where 


n\ 


plk  =  o  dmPi  (m)  (25) 

n\- 1 

Computational  efficiency  like  that  for  the  bin-to-bin  cross  sections  is  achieved  for  the 

(21  +  1) 

Legendre  moments  by  pre-computing  and  storing  an  array  of - Pf  k  for 


/  =  OK  L  and  k  =  IKK,  for  a  choice  of  L  and  K. 


II.4.3:  Multi-Group  Monte  Carlo  Transport 

In  addition  to  the  fully  discretized  transport  methods,  multi-group  Monte  Carlo 
transport  could  also  use  the  PAX  cross  sections  by  creating  a  tabular  cumulative 
distribution  function  to  invert  by  table  search  and  interpolation.  The  conventional 
method  for  generating  group-to-group  cross  sections  for  multi-group  Monte  Carlo 
creates  a  strictly  positive  representation  of  the  scatter  cross  section  from  a  truncated 
Legendre  expansion  using  a  maximum  entropy  method  (5,  1 1).  Sixteen  equally-likely 
intervals  of  m  are  then  created  from  this  strictly  positive  representation.  As  shown  in 
Figure  1,  however,  the  group-to-group  cross  sections  can  have  several  separated 
regions  of  zero  value.  Neither  the  16  equal-likelihood  intervals  nor  the  strictly  positive 
maximum  entropy  method  accurately  approximates  these  zero-value  regions. 

Therefore,  I  expect  the  PAX  cross  sections  would  be  more  accurate  than  the 
conventional  method  for  creating  group-to-group  Monte  Carlo  cross  sections  (although 
this  research  is  left  for  a  future  effort). 

1 1. 5:  Failure  of  Point- Wise  Evaluation  of  the  Group-to-Group  Cross 
Sections 

The  PAX  cross  sections  are  a  finite  volume  approach  rather  than  a  numerical 
quadrature  based  on  interpolation  between  point  values,  such  as  composite  midpoint: 
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An  attempt  was  made  to  evaluate  the  group-to-group  cross  sections  at  a  grid  of 
points  for  use  within  a  numerical  quadrature  based  on  interpolation.  Point-wise 

evaluation  of  the  group-to-group  cross  section  fails  due  to  poles  in  s  (rt^  )  for  the 

level  inelastic  scatter  mechanism  (where  is  the  cosine  of  the  scatter  angle  in  the 


laboratory  frame  of  reference).  The  level  inelastic  scatter  mechanism  can  have  an 
ENDF/B-VI  representation  in  the  center  of  mass  frame  of  reference.  Evaluating  the 
cross  section  in  the  laboratory  frame  (the  frame  of  the  transport  problem)  involves 
transforming  the  distribution  function  for  the  cosine  of  the  scatter  angle  into  the 
laboratory  frame  using 
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where  f  represents  the  distribution  function  for  the  special  case  of  the  angular 
distribution  described  in  section  III.  1.4,  and  mCM  is  the  cosine  of  the  scattering  angle  in 

the  center  of  mass  frame  of  reference.  For  the  level  inelastic  scattering  mechanism,  this 
transformation  is  infinite  as  approaches  1  and  the  incident  energy  of  the  neutron 


approaches  the  energy  deficit.  The  result  is  that  the  incident  energy  integration  in 
equation  (26)  fails  to  converge. 

Numerically  evaluating  the  PAX  cross  sections 
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is  not  a  problem  because  s  (E  E,nfM  ^  is  well-behaved,  and 


dm. 


CM 


dnij . 


appear  because  H  is  a  point  function,  not  a  distribution  function. 


does  not 
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III.  Implementation 

The  numerical  approximations  to  the  group-to-group  scatter  cross  sections 
using  the  scatter  cross  section  operator  defined  in  chapter  2  have  been  implemented  in  a 
computer  code.  This  computer  code  uses  ENDF/B-VI  data  to  generate  the  three 
different  numerical  approximations  using  the  scatter  cross  section  operator.  Although  a 
computer  code  could  have  been  written  to  implement  only  the  new  PAX  cross  sections, 

I  decided  to  implement  a  single  code  that  could  perform  all  three  different  group-to- 
group  cross  section  approximations:  Legendre  moments,  bin-to-bin  cross  sections,  and 
PAX  cross  sections,  because  the  new  PAX  cross  sections  would  have  to  be  validated 
against  existing  codes  that  only  output  either  bin-to-bin  cross  sections  or  Legendre 
moments. 

My  code  philosophy  balances  efficiency,  robustness,  and  code  readability.  The 
programming  language  FORTRAN  90/ 95  has  self-documenting  features  such  as 
extended  variable  names  and  modular  structure.  Coding  decisions  were  made  to 
emphasize  readability  for  the  purpose  of  debugging.  Therefore,  the  improvements  in 
efficiency  in  the  computation  are  a  consequence  of  the  algorithms  discussed  in  this 
chapter  and  the  cross  section  operator  discussed  in  chapter  2.  Code  optimization  is  left 
to  the  compiler. 

This  chapter  presents  several  algorithms  used  in  the  new  computer  code.  These 
algorithms  use  some  FORTRAN  90/ 95  key  words.  These  key  words  are  summarized 
in  appendix  B,  and  are  printed  in  an  italic,  mono-spaced  font,  e.g.  El  self. 

ENDF/B-VI  data  is  used  directly  so  that  no  intermediate  approximations  are 
used  prior  to  calculating  the  group-to-group  cross  sections.  DelGrande  and  Mathews 
(7)  used  an  intermediate  output  from  the  NJOY  code  called  A  Compact  ENDF  (ACE) 
file  (13).  The  NJOY  code  uses  approximations  to  generate  the  ACE  files.  These 
approximations  introduce  errors  that  cannot  be  controlled  by  the  user  of  the  new 
computer  code.  NJOY  has  its  own  residual  constraints  from  earlier  assumptions  about 
computational  power  and  resources  because  it  was  developed  during  the  1970s. 

III.  i :  Cases  for  Secondary  Energy  and  Angular  Distributions 

In  general,  ENDF/B-VI  represents  scatter  cross  sections  as 


20 


(Em  E,m)=  s  (Et)„(Et)fjoM  (m,E\Ef), 


(30) 


where  the  joint  distribution  is  normalized  as 

odEod mf joun  fa  E  \E  & =  (3 1 ) 

and  the  multiplicity  function,  n(E  ty,  represents  the  number  of  secondary  neutrons 

created  as  a  function  of  the  incident  energy.  The  multiplicity  function  is  often  used  for 
fissionable  isotopes  but  can  also  be  used  to  represent  the  number  of  secondary  neutrons 
in  reactions  such  as  (n,  2n).  With  the  inclusion  of  the  multiplicity  function  in  the  scatter 
cross  section  by  ENDF/B-VI,  the  fission  cross  section  is  a  mechanism  to  be  used  within 
the  scatter  cross  section  operator. 

The  joint  distribution  can  be  expressed  in  three  ways, 

j  f(m\Et)g(E\Et) 

fjoint  faE  \E$)=  jg(£  \Efy(m\E<tiE).  (32) 

| /  (m\E  fyg(E  \E  i  m) 

Each  of  these  three  cases,  have  approximations  of  f  and  g,  and  those  approximations  are 
called  laws  by  ENDF/B-VI  notation.  Elastic  and  level-inelastic  scatter  are  treated 
together  as  a  special  case. 

All  three  of  these  cases  can  be  represented  in  ENDF/B-VI  in  either  the  center  of 
mass  reference  frame  or  the  laboratory  reference  frame.  If  the  ENDF/B-VI  data  is 
represented  in  the  center  of  mass  reference  frame,  then  the  scatter  cross  section 
operator  is  transformed  from  the  laboratory  reference  frame  to  the  center  of  mass 
reference  frame.  Otherwise,  if  the  ENDF/B-VI  data  is  represented  in  the  laboratory 
reference  frame,  then  the  scatter  cross  section  operator  is  not  transformed. 

The  cross  section  function,  s  (E  ty,  may  be  quite  smooth  over  large  energy 

ranges  or  it  can  vary  rapidly  in  regions  where  resonances  are  present  for  a  given 
material.  The  rapidly  varying  resonance  regions  are  represented  using  parameterized 
functions  with  tabulated  parameters.  Other  regions  are  represented  by  interpolation 
functions  using  tabulated  values. 

The  scatter  cross  section  operator  is  restated  for  convenience  with  the  above 
approximation  of  the  continuous  scatter  cross  section  from  ENDF/B-VI  as 


21 


(33) 


g= 


j?  asi  $r 

5  dE  05  (E  fyi  (E  0)F  (E  0)|  5  dE  |q  dmfjoint  (™>E  \E  0)g|± 

DEgt _ BD^g  *  1 _ 5 

5  dE  fiF  (E 


D  E 


gt 


The  explicit  nesting  of  the  integrals  in  the  scatter  cross  section  operator  is  used  to  show 
the  different  dependencies  of  the  cases. 

Three  cases  for  the  joint  distribution  are  presented.  Then,  elastic  and  level 
inelastic  scatter  are  presented  as  an  important  special  case.  The  mechanisms  that 
typically  use  each  of  the  cases  are  presented.  The  actual  ENDF/B-VI  laws  are  stated  in 
appendix  F. 


III.  1 . 1 :  Separable  Energy  and  Angular  Distributions 


The  assumption  in  the  first  case  is  the  separability  of  the  secondary  energy 
distribution  and  the  scatter  angle  distribution.  It  is 

fjouu  (”* E  /('T  \E*>  (»*) 


This  distribution  is  often  used  for  representations  such  as  fission  or  (n,  2n) 
reactions.  This  representation  is  only  chosen  when  the  reactions  produce  at  least  three 
secondary  particles,  ie.  two  neutrons  and  a  remaining  nucleus.  This  three-body  problem 
results  in  a  weak  dependence  between  the  angular  distribution  and  the  secondary 
energy  distribution  of  the  secondary  particles.  The  ENDF/B-VI  evaluators  can  then 
reasonably  separate  the  joint  distribution  into  two  separate  distributions. 

With  separable  energy  and  angular  distributions,  the  integrals  in  the  scatter 
cross  section  operator  can  be  rearranged  as 
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Calculating  a  mechanism  with  a  separable  energy  and  angular  distribution  is  the 
most  efficient  case  because  the  integration  requires  only  two-level  nesting  of 


22 


quadratures.  Section  III.4.1  presents  the  algorithm  I  use  to  calculate  the  two  innermost 
integrations. 


III.  1.2:  Angular  Distribution  Dependent  on  Secondary  Energy 

The  second  case  has  an  angular  distribution  dependent  on  the  secondary  energy. 
This  case  is 


4,„,  (»*£  S(E  \E^f('n\EiE).  (36) 


This  case  can  be  used  for  any  (n,  n+product)  reactions,  where  the  product  can  be 
a  neutron,  multiple  neutrons,  a  proton,  an  alpha  particle,  etc.  It  is  used  when  the 
secondary  energy  and  the  angular  distribution  are  more  strongly  correlated.  It  can  also 
be  used  when  the  ENDF/B-VI  evaluators  have  so  little  information  that  only  a  rough 
assumption  about  the  secondary  energy  and  angular  distributions  can  be  made.  The 
rough  approximation  to  the  angular  and  energy  correlation  may  be  due  to  the  large 
number  of  products  produced  during  the  reaction.  For  both  possibilities,  the  second 
case  can  represent  the  appropriate  ENDF/B-VI  laws.  The  integrals  in  the  scatter  cross 
section  operator  can  be  rearranged  as 
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This  case  takes  the  most  computational  time  due  to  the  three-level  nesting  of  the 
quadrature. 


III.  1.3:  Secondary  Energy  Distribution  Dependent  on  Scatter  Angle 

The  last  case  is  used  for  only  one  law  in  ENDF/B-VI:  the  laboratory  angle- 
energy  law.  This  case  is 

f joint  (n%E  \E  ^)=  / (m\E  t )g(E  \E  £m)-  (38) 

In  rare  circumstances,  the  ENDF/B-VI  evaluators  can  only  accurately 
determine  the  secondary  energy  of  the  products  of  the  reaction  as  a  function  of  angle,  as 
in  a  laboratory  experiment.  Then,  the  secondary  energy  distribution  is  clearly 
dependent  on  the  scatter  angle  and  this  appropriate  law  is  chosen. 
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The  use  of  the  laboratory  angle-energy  law  in  ENDF/B-VI  is  so  rare  that  I  did 
not  find  an  example  that  used  it.  Nevertheless,  my  code  does  support  this  option  by 
rearranging  the  integrals  in  the  scatter  cross  section  operator  as 
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Presumably,  this  is  also  computationally  expensive  because  of  the  three-level 
nesting  of  the  quadrature. 


III.  1.4:  Elastic  and  Level  Inelastic  Scatter  Cases 

The  special  case  for  the  third  case  is  elastic  and  level  inelastic  scatter 
mechanisms.  In  these  mechanisms,  the  secondary  energy  of  the  neutron  is  uniquely 
determined  by  the  incident  energy  and  the  angle  of  scatter,  with  the  formula 

fJoint  (rn,E  | Efy=  f(m\E$)d(E  -  Es  (E$mfr  (40) 


where  E,  is  determined  by  conservation  of  energy  and  momentum  as 
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and  where  is  “the  excess  of  kinetic  energy  of  the  product  particles  over  that  of  the 
original  particles”  (l  l),  and  A  is  the  ratio  of  the  mass  of  the  target  nucleus  to  the  mass  of 
the  neutron. 

The  integrals  in  the  scatter  cross  section  operator  can  be  rearranged  as 
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where  the  secondary  energy  group  boundaries,  Eg  and  Eg~i,  are  included  in  the  bounds  of 
the  integration  with  respect  to  m.  The  bounds  can  be  obtained  from  equation  (41) 


given  E  ft,  E,  and  0. 
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The  special  case  of  elastic  and  level  inelastic  scatter  is  the  most  frequently  used 
and  typically  contributes  the  most  to  the  group-to-group  scatter  cross  sections. 
Although  finding  the  limits  of  integration  of  the  cosine  of  the  scattering  angle  takes 
extra  computational  time,  the  overall  computational  effort  is  comparable  to  the  first  case 
(section  III.  1.1 )  of  separable  energy  and  angular  distributions.  Additionally,  because  the 
neutron  can  always  scatter  elastically  at  any  incident  energy,  this  special  case  is  always 
used  for  each  incident  energy  group. 

1 1 1. 2:  Quadratures 

Two  adaptive  quadratures,  Simpson  and  Gauss-Simpson,  were  used  to  perform 
the  needed  nested  integrations.  Both  adaptive  quadratures  used  FORTRAN  90/ 95 
recursion  to  subdivide  intervals  as  necessary  to  achieve  user-set  error  tolerances.  Two 
other  quadrature  methods  were  also  considered,  but  were  discarded — Romberg 
integration  (4)  and  IMSL  (5)  adaptive  integration  routines. 

III. 2.1:  Characteristics  of  the  Integrand 

Section  III.  1  presented  four  cases  of  nested  integrations.  Each  of  the  nested 
integrations  uses  data  from  ENDF/B-VI.  The  data  from  ENDF/B-VI  is  associated 
with  laws  that  describe  the  use  of  either  tabulated  data  or  tabulated  parameters  for 
reconstruction  of  the  data.  For  either  type  of  ENDF/B-VI  data,  the  characteristics  of 
the  integrands  in  section  III.  1  are  similar: 

1.  Finite  number  of  discontinuities  in  the  integrand 

2.  Finite  number  of  discontinuities  in  the  first  derivative  of  the  integrand 

3.  Localized,  non-polynomial  behavior 

4.  Finely  partitioned  domain  to  handle  the  characteristics  in  #1  and  #2  above 

A  quadrature  with  a  robust,  adequate,  and  practical  implementation  was  needed 
to  address  the  characteristics  of  the  integrand.  The  discontinuities  in  the  integrand  and 
the  first  derivative  of  the  integrand  require  a  relatively  fine  initial  mesh  for  any  of  the 
integrations  described  in  section  III.  1 .  With  a  relatively  fine  initial  mesh,  a  modest 
order  numerical  quadrature  is  sufficient.  The  localized,  non-polynomial  behavior 
requires  an  adaptive  method.  An  open  quadrature  method  does  not  evaluate  the 
integrand  at  the  endpoints  of  the  domain,  but  a  closed  quadrature  method  uses  at  least 
the  endpoints  of  the  domain  of  integration.  And,  to  avoid  data  and  code  complexity  to 
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use  only  a  closed  method,  both  open  and  closed  methods  were  used  when  appropriate. 
An  algorithm  of  modest  order,  adaptive,  and  with  either  open  or  closed  methods 
addressed  each  of  the  characteristics  of  the  integrand. 


III. 2.2:  Romberg  Integration 

The  first  numerical  quadrature  discarded  was  a  Romberg  automatic  integration, 
given  in  numerical  methods  texts  such  as  Burden  and  Faires  (4).  The  Romberg 
numerical  method  was  discarded  for  two  reasons,  global  error  testing  and  closed  end¬ 
point  method.  The  Romberg  method  uses  a  global  convergence  test.  I  used  the  relative 
error 
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If  the  entire  integration  did  not  pass  the  relative  error  tolerance,  then  another  Romberg 
iteration  was  performed.  Because  each  Romberg  iteration  is  twice  as  expensive  as  the 
previous  iteration,  this  method  was  slower  to  converge  than  the  adaptive  quadrature 
methods  used.  In  practice,  the  global  error  in  the  integration  could  be  dominated  by 
only  one  or  two  sub-domain  pieces.  Despite  this,  Romberg  subdivides  all  the 
subintervals. 

Romberg  integration  is  a  succession  of  composite  trapezoid  quadratures  with  an 
Aitken  extrapolation  for  the  next  approximation  to  the  integral.  The  trapezoid 
numerical  quadrature  is  a  closed  quadrature,  which  is  not  useful  for  several  of  the  nested 
quadratures  required  in  evaluating  the  integrals  in  section  III.  1 .  One  could  use  a 
Romberg  scheme  based  on  composite  midpoint,  which  is  an  open  rule,  but  it  is  even  less 
efficient.  In  order  to  reuse  function  evaluations,  it  requires  3n  evaluations  for  n  levels  of 
subdivision;  composite  trapezoid  requires  only  2n. 


III. 2. 3:  IMSL  Adaptive  Integration  Routines 

IMSL’s  adaptive  quadratures  (5)  were  not  used  for  two  reasons.  The  first  reason 
is  that  they  use  Gauss-Kronrod  quadrature  in  each  subinterval;  thus  the  number  of 
points  used  was  excessive — 17  points  in  each  subinterval.  Usually,  the  partitioning  of 
the  domain  to  account  for  discontinuities  was  sufficient  for  the  relative  error  tolerance 
to  be  passed  with  a  far  more  modest  number  of  points,  such  as  those  used  in  either 
Simpson  or  Gauss-Simpson  adaptive  methods. 
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The  second  reason  that  the  IMSL  routines  were  not  used  is  the  many  required 
different  subroutines  and  variations.  The  IMSL  integration  routines  are  written  in 
FORTRAN  77  and  pre-compiled,  and  consequently  are  not  recursive.  Without 
recursion,  the  same  IMSL  library  routine  cannot  be  used  for  nested  integrals. 
Therefore,  the  IMSL  routines  were  discarded  due  to  a  lack  of  flexibility  and  excessive 
number  of  quadrature  points. 

III. 2.4:  Adaptive  Simpson  Integration 

The  limits  of  integration  and  the  integrand  are  assumed  to  be  separate 
subroutines  to  either  adaptive  integration  quadrature.  Given  the  limits  of  integration, 
any  places  in  the  domain  where  the  integrand  is  discontinuous  or  its  first  derivative  is 
discontinuous  must  be  found  and  tabulated.  The  initial  tabulation  is  used  as  a  mesh  for 
the  adaptive  quadratures.  The  adaptive  quadratures  then  integrate  each  of  the  mesh 
pieces  individually.  Each  of  the  adaptive  methods  performs  two  quadratures  on  the 
subinterval.  If  the  two  quadratures  pass  a  relative  error  test  for  convergence,  then  the 
integration  on  the  subinterval  is  complete;  otherwise,  the  subinterval  is  subdivided  and 
each  sub-piece  is  integrated  and  tested  for  convergence  individually. 

The  adaptive  Simpson  routine  (Algorithm  l)  is  logically  correct.  In  the  code,  it 
is  implemented  with  additional  arguments  in  order  to  avoid  redundant  function 
evaluations  and  quadrature  calculations.  The  algorithm  for  the  adaptive  Simpson 
method  is  covered  first  because  the  adaptive  Simpson  method  is  embedded  in  the 
adaptive  Gauss-Simpson  method. 

Algorithm  1:  Adaptive  Simpson 

Input:  Xq,  Xj 

Output:  Q 

Dx  =  Xj  -  xQ 

Q0  =  (f(xo)+  4/(x  0  +  Dx/2)  +  /(*  l))Dx/3 

Ql  =  (/(*o)+  4/(x0  +  Dx/4  )+/(xQ  +  l/2Dx))Dx/6 

Qr  =  (/(*  o  +  DV2)+  4/(xo  +  3Dx/4)+  /(xt))Dx/6 

Q\  =  Qr  +  Ql 

Q0 1  £  erd  +  Qq)  )  Then 
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0=0! 

Else 

X0  +  X, 

Left  integral  =  Adaptive  Simpson  (xQ,  xQ  +  - - - ) 

x0  +  x. 

Right  integral  =  Adaptive  Simpson  (x()  +  - — - ,  Xj ) 

Q  —  Left  integral  +  right  integral 
End  If 


III. 2. 5:  Adaptive  Gauss-Simpson  Integration 

The  adaptive  Gauss-Simpson  initially  uses  a  two  point  Gauss-Legendre 
quadrature.  If  the  relative  error  tolerance  is  not  passed,  then  the  Gauss  subroutine  calls 
two  different  Gauss-Simpson  subroutines  for  each  portion  of  the  integrand. 


Algorithm  2:  Adaptive  Gauss 
Input:  x0,  Xj 
Output:  Q 
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03® 


*-* 

2 


x 


X„  +  X  .  , 
_  0  m  id 


m  id, 


X 


X  .  ,  +  X, 
_  mid  1 


mid, 

Qr  =  h 


2  2 
as  as 


m  id. 


h  6  as  h  & 

j  &  mid  +  /V# 


03  £> 


03® 


h  6  as  h 

03 ' 7  §*"‘*2  +  03  S 


0!  -  0*  + 

If  ( 2^!  -  0O|  £  eg/  (0!  +  0O)  ;  Then 

0=0! 

Else 
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Xq  +  X. 

Left  integral  =  Adaptive  Left  Gauss-Simpson  (xQ,  xQ  +  - — - ,QL  ) 

x0  +  x. 

Right  integral  =  Adaptive  Right  Gauss-Simpson  (xQ  +  - — - ,  Xy,QR  ) 

Q  —  Left  integral  +  right  integral 
End  If 


Algorithm  3:  Adaptive  Left  Gauss-Simpson 
Input:  Xq,  X y,  Qprev 


Output:  Q 

X,  -  Xr 

h 


1  ^0 


mid 


mid. 


XQ  +  Xy 


Xn  +  X 
u  m  id 


Qr  =  h 


as  as 


m  id. 


h  6  as  h  c& 


V3^ 


V3® 


Dx  =  x,  -  x  . , 

1  mid 


Qr  =  (/'(Ym,/)+  4f(xmid  +  Dx/2)+/(x1))Dx/3 


0!  - 

If  (  2\Qr  Qo\ £  erfi/ (& +  fio) ; Then 
Q  =  QX 

Else 

Xq  +  X, 

Left  integral  =  Adaptive  Left  Gauss-Simpson  (xQ,  xQ  +  - — - ,QL  ) 

x0  +  x. 

Right  integral  =  Adaptive  Simpson  (xQ  +  - — - ,  Xy ) 

Q  =  Left  integral  +  right  integral 
End  If 


Algorithm  4:  Adaptive  Right  Gauss-Simpson 
Input:  XQ,  Xy,  Qprev 
Output:  Q 

Xn  +  X, 

X  . ,  =  - - 

mid  2 

Dx  =  x  . ,  -  xn 

mid  U 

=  (/(xo)+  4/(xi  +  D*/2)+  fix  mi  j)°x/3 

/2  =  Xl  '  X° 

4 
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X  .  ,  +  X, 
mid _ 1 


mid, 

Qr  =  h 


2  2 
as  as 


m  id~ 


h  i+  & 


+ 


h  <M 


03  S  "  fmidi  '  03# 


Qr  +  Ql 

If  (  2 10!  -  0O |  £  e.e/  (0!  +  00 )  ;  Then 

0=0! 

Else 

xo  + 

Left  integral  =  Adaptive  Simpson  (xQ,  +  - - — -) 

x0  +  x. 

Right  integral  =  Adaptive  Right  Gauss-Simpson  (xQ  +  - - - ,  Xp0^  ) 

Q  —  Left  integral  +  right  integral 
End  If 


The  adaptive  Simpson  integration  is  more  efficient  than  the  adaptive  Gauss- 
Simpson  because  all  of  the  points  in  the  integration  can  be  reused  if  additional 
subdivision  is  needed.  This  efficiency  can  improve  the  runtime  of  the  calculation  by 
approximately  a  factor  of  two. 

The  adaptive  Gauss-Simpson  integration  method  is  used  because  it  is  an  open 
rule.  A  two-point  Gauss-Legendre  quadrature  method  was  used  with  composite 

Simpson  because  both  methods  have  error  (9(Dx)4  .  Having  more  points  in  the  Gauss- 
Legendre  method  would  decrease  efficiency  because  the  points  are  discarded  upon 
recursion. 

III. 3:  Main  Program  and  Incident  Energy  Integration 

The  portions  of  the  algorithm  that  do  not  involve  the  scatter  cross  section 
operator  reside  in  the  main  program.  These  include  file  input/ output,  looping  through 
the  incident  energy  groups,  looping  through  the  different  mechanisms,  and  looping 
through  the  secondary  energy  groups.  In  the  loop  for  the  secondary  energy  group,  a 
subroutine  is  called  to  perform  the  outer-most  integral,  the  incident  energy  integral,  for 
the  scatter  cross  section  operator  for  all  four  cases  presented  in  section  III.  1 

Algorithm  5:  Main  program 

Call  GetENDFData 

Do  g'  =  1,  G  (number  of  energy  groups) 
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Do  mech  =  1,  M  (number  of  mechanisms) 

If  (  mechanism  does  not  occur  )  Cycle 
Call  GetLimits(  giow,  ghigh  ) 

DO  g  =  glow,  ghigh 

Call  IncidentEnergy(  integral  value(g,  mech) ) 
End  Do 

Accumulate  the  sum  of  the  integrals 
End  Do 

Call  OutputIntegral(g'-to-g  integral ) 

End  Do 


Because  the  incident  energy  integration  has  no  discontinuities  in  the  integrand, 
it  uses  the  recursive,  adaptive  Simpson  method  detailed  in  section  III. 2  for  its  efficiency. 
The  incident  energy  mesh  is  constructed  with  special  care  because  of  the  many  places 
within  the  limits  of  integration  for  the  incident  energy  group  that  its  integrand,  has,  or 
could  have,  discontinuities  in  the  first  derivative.  These  discontinuities  arise  from  the 
following  places: 

1.  User-specified  boundaries  in  the  energy  group  structure 

2.  Tabulation  of  s  (E  from  ENDF 

3.  Tabulation  of  n  (E  from  ENDF 

4.  Tabulation  of  F  (E  from  user-set  parameters 

5.  Tabulation  of  / (in\E  0)  or  / (m\E  %E )  from  ENDF 

6.  Discontinuous  tabulation  of  either  PAX  cross  sections  or  direct  calculation  of 

bin-to-bin  cross  sections  with  h  ( m )  when  using  elastic  or  level-inelastic 

scatter  mechanisms  because  of  implicit  secondary  energy  group  dependence 

7.  Tabulation  of  either  g^E  \E  or  g^E  \E  m ^  from  ENDF 

Given  an  incident  energy  mesh  with  mesh  points  at  all  of  the  discontinuities  in 
the  first  derivative,  the  algorithm  for  the  incident  energy  mesh  calls  the  recursive, 
adaptive  Simpson  method  given  in  Algorithm  1  above  for  each  mesh  interval. 
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III. 4:  Joint  Angular  and  Secondary  Energy  Distribution  Algorithms 

The  joint  angular  and  secondary  energy  distributions  use  the  three  cases  and  the 
special  case  given  in  section  III.  1 .  The  four  cases  use  slightly  different  algorithms  and  I 
present  them  individually. 

III.4.1:  Separable  Angular  and  Secondary  Energy  Distribution 

The  first  case  that  is  detailed  in  section  III.  1.1,  separable  angular  and  secondary 
energy  distributions,  uses  two  separate  (as  opposed  to  nested)  integrations.  The 
separable  energy  integration  is  performed  using  closed-form  solutions  for  integrals  of 
the  functions  using  the  tabulated  parameters  from  ENDF/B-VI.  The  separable  angular 
distribution  is  performed  using  the  adaptive  Gauss-Simpson  method  presented  in 
section  III. 2  above.  The  PAX  cross  section,  equation  (20),  has  discontinuities  in  the 
cosine  of  the  scatter  angle  mesh  from  the  two  Heaviside  functions.  Because  I  chose  to 
include  all  of  the  different  approximations  for  generating  the  group-to-group  scatter 
cross  sections  from  the  cross  section  operator  in  one  computer  code,  these 
discontinuities  require  use  of  the  adaptive  Gauss-Simpson  method. 

III. 4. 2:  Angular  Distribution  Dependent  on  Secondary  Energy 

The  second  case  presented  in  section  III.  1.2  has  an  angular  distribution 
dependent  on  the  secondary  energy.  This  algorithm  uses  two  nested,  adaptive  Gauss- 
Simpson  calls.  The  secondary  energy  integration  can  have  discontinuities  in  the 
tabulation  for  the  secondary  energy.  Therefore,  the  open  adaptive  Gauss-Simpson 
method  was  an  obvious  choice.  The  integration  for  the  angular  distribution  used  the 
Gauss-Simpson  method  for  the  same  reason  as  given  in  the  separable  case. 

III. 4.3:  Secondary  Energy  Distribution  Dependent  on  Scatter  Angle 

The  third  case  presented  in  section  III.  1.3,  in  which  the  secondary  energy 
distribution  is  dependent  on  the  scatter  angle,  uses  different  integration  methods  for  the 
outer  and  inner  integrals.  The  outer  integral,  for  the  angular  distribution,  uses  an 
adaptive  Gauss-Simpson  method  as  given  in  the  separable  case.  The  inner  integral  uses 
closed-form  solutions  for  the  integrations  of  the  interpolating  functions  using  the 
tabulated  values  of  the  secondary  energy  given  by  the  ENDF/B-VI,  laboratory  angle- 
energy  law. 
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III. 4.4:  Elastic  and  Level  Inelastic  Scatter 


The  elastic  and  level  inelastic  scatter  special  case  described  in  section  III.  1.4  uses 
the  adaptive  Gauss-Simpson  method  for  the  same  reason  given  in  section  III.4.1  for  the 
separable  angular  distribution.  There  is  no  secondary  energy  group  integration  given 
in  equation  (42).  To  account  for  the  secondary  energy  group,  a  subroutine  to  calculate 
the  cosine  of  the  scatter  angle  limits  of  integration  dependent  on  the  secondary  energy 
group  is  added  prior  to  the  adaptive  Gauss-Simpson  integration. 

To  handle  the  laboratory  frame  transformation  to  the  center  of  mass  frame,  an 
additional  subroutine  was  required  to  transform  the  boundaries  of  the  mesh  used  for  the 
cosine  of  the  scatter  angle  for  either  the  PAX  cross  sections  in  equation  (20)  or  for 
direct  evaluation  of  the  bin-to-bin  cross  sections  in  equation  (17).  A  point-wise 
transformation  was  also  required  for  the  evaluation  of  the  Legendre  moments  in 
equation  (19).  The  same  equation  used  for  the  transformation  can  be  derived  from 
equation  (41)  and  is  described  in  detail  in  the  NJOY  documentation  (13). 

1 1 1. 5:  Data  Structure 

ENDF/B-VI  allows  cross-section  evaluators  to  store  data  in  any  of  numerous 
different  ways.  Compromises  between  memory  requirements  and  efficiency  in  accessing 
the  data  within  my  code  required  some  ingenuity  and  the  use  of  FORTRAN  90/ 95 
allocatable  arrays  for  dynamic  memory  usage. 

These  arrays  are  allocated  at  runtime  to  store  the  largest  extent  of  each  of  the 
dimensions  needed  for  a  mechanism’s  distribution.  For  example  10B  has  several  level 
inelastic  scatter  mechanisms  available  as  well  as  an  elastic  scatter  mechanism.  Each  of 
the  mechanisms  has  potentially  different  data  storage  and  representation.  Therefore,  I 
used  derived  types  to  allow  for  all  of  the  available  choices  for  the  angular 
distributions.  This  derived  type  contains  the  integer  flags  necessary  to  identify 
the  appropriate  representation  (one  of  the  ENDF/B-VI  laws  listed  in  appendix  F). 

Then,  an  array  is  allocated  for  each  angular  distribution  representation  to  the  largest 
extent  required  by  any  of  the  mechanisms.  Finally,  the  data  is  read  from  the  input 
ENDF/B-VI  fde  and  fdled  into  the  appropriate  representation  array.  One  of  the 
consequences  of  this  approach  is  that  the  code  must  allocate  extra  integer  arrays  that 
contain  the  boundaries  for  each  of  the  data  arrays.  To  demonstrate  how  the  arrays  are 
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allocated  and  how  the  data  is  read  into  these  arrays,  a  sample  of  a  pseudo-code  for  one 
type  of  storage  is  provided  in  Algorithm  6. 


Algorithm  6:  Pseudo-code  for  reading  and  storing  ENDF/B-VI  data 

Open  (  ENDF  file  ) 

Do 

Read  (  Character  String,  line  number  ) 

If  ( line  number  ==  exit  line  number  )  Exi  t 
End  Do 
Do 

Read  (  mechanism  number,  ENDF  law  ) 

If  (  mechanism  number  /  =  valid  mechanism  number  )  Cycle 
Select  Case  (  ENDF  law  ) 

Case  (  Legendre  expansion  ) 

Read  (  number  of  energy  mesh,  number  of  Legendre  moments  ) 
max  number  energy  mesh  =  Max  (  new  number,  previous  max  ) 
max  number  moments  =  Max  (  new  number,  previous  max  ) 

End  Case 

Read  (  next  section  number  ) 

Iff  next  section  number  /  =  this  section  number  )  Exi  t 
End  Do 

Close  (  ENDF  file  ) 

Allocate  (  energy  data(max  incident  mesh)  ) 

Allocate  (  moment  data(0:max  number  moments,  max  incident  mesh)  ) 
Open  (  ENDF  file  ) 

Do 

Read  (  mechanism  number,  ENDF  law  ) 

If  (  mechanism  number  /  =  valid  mechanism  number  )  Cycle 
Select  Case  (  ENDF  law  ) 

Case  (  Legendre  expansion  ) 

Read  (  number  of  energy  mesh,  number  of  Legendre  moments  ) 
Do  i  =  1,  number  of  energy  mesh 
Read  (  energy  data(i)  ) 

Read  (  moment  data(0:number  of  Legendre  moments)  ) 

End  Do 
End  Case 

Read  (  next  section  number  ) 

Iff  next  section  number  /  =  this  section  number  )  Exi  t 
End  Do 


34 


IV.  Cross  Section  Code  Validation 


I  validated  the  cross  section  code  by  comparing  its  output  with  the  results  of  two 
other  codes:  NJOY  and  DelGrande’s  Monte  Carlo  discrete  elements  cross  section 
codes.  Several  online  ENDF/B-VI  plotting  sites  were  also  used  for  initial  validation 
checks  (19,  21).  NJOY  was  used  to  validate  the  accurate  evaluation  of  the  secondary 
energy  and  angular  distributions  as  well  as  the  nested  integrations.  The  online  sites 
were  used  to  directly  validate  the  values  produced  by  my  code  of  the  cross  section  data 
recreated  from  ENDF/B-VI.  Comparing  the  new  computer  code  to  DelGrande’s  Monte 
Carlo  discrete  elements  code  validated  the  discrete  elements  approximation. 

IV.  i:  Direct  Validation 

Validation  began  by  confirming  that  the  input  data  had  been  read  correctly  from 
the  ENDF/B-VI  file.  After  checking  the  input,  it  was  important  to  verify  that  the  data 
was  being  faithfully  represented  and  recreated.  When  dealing  with  the  many  different 
interpolation  functions,  ENDF/B-VI  laws,  and  different  reference  frames,  the  faithful 
reconstruction  of  the  cross  section  curves  was  essential. 

Direct  validation  of  the  cross  section  curves  involved  examining  and  comparing 
many  different  plots  for  many  different  isotopes  and  mechanisms.  The  plots  were  then 
superimposed  to  determine  whether  or  not  the  curves  generated  by  my  program 
matched  those  at  respected  online  sites  such  as  the  ENDF/B-VI  site  (19)  or  the  T2- 
Nuclear  site  (21),  which  is  a  site  run  by  the  group  within  the  Los  Alamos  national 
laboratory  that  maintains  and  distributes  the  NJOY  code. 

All  of  the  different  resonance  region  parameterizations  with  Doppler  broadening 
of  the  cross  section  and  the  tabulated  cross  sections  were  checked.  The  joint 
angular/ secondary  energy  distributions  were  left  to  the  NJOY  validation  portion. 
Figure  7  shows  an  example  of  one  of  the  validations  performed  using  the  ENDF/B-VI 
site  (19).  The  curve  is  the  elastic  scatter  cross  section  of  56Fe  using  the  Reich-Moore 
resonance  parameterization,  which  is  one  of  the  most  complicated  cross  section  curves 
represented  in  ENDF/B-VI  files.  To  the  limit  of  the  benchmark  data,  this  calculation  is 
correct. 
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Incident  neutron  energy  (eV) 

Figure  7:  Comparison  of  download  data  to  reconstructed  data 

S6Fe  elastic  scatter  cross  section  reconstruction  of  Reich-Moore 
resonance  region  parameterization  broadened  to  300K 

IV. 2:  NJOY  Validation 

Although  the  direct  validation  of  the  cross  section  curves  is  a  powerful  tool,  it 
cannot  validate  everything  in  my  new  computer  code.  The  validation  of  the 
reconstruction  and  integration  of  the  angular  and  secondary  energy  distributions,  as 
well  as  the  integration  of  the  cross  section,  used  the  NJOY  code.  Two  methods  using 
NJOY  were  employed — the  development  environment  that  contains  a  powerful 
debugging  tool,  and  the  output  of  the  NJOY  calculation. 

IV. 2.1:  NJOY  Validation  Using  the  Development  Environment 

The  Compaq  Visual  FORTRAN  development  environment  (5)  lias  two 
debugging  features  that  I  used  for  validation.  Both  debugging  tools  are  used  while 
performing  step  debugging,  a  runtime  environment  that  allows  a  programmer  to  step 
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through  the  execution  of  a  computer  code  line  by  line.  The  first  feature  allows  the 
programmer  to  view  the  value  of  any  variable  while  step  debugging.  The  second 
feature  allows  the  programmer  to  change  the  value  of  any  variable  during  step 
debugging.  Using  these  two  tools,  it  was  possible  to  enter  a  specific  value  for  any  of  the 
different  distributions  in  the  new  computer  code.  The  calculated  result  of  using  the 
entered  value  could  then  be  compared  to  the  NJOY  calculated  result  while  both 
computer  codes  were  running  in  a  step  debugging  mode.  Not  only  did  this  greatly  aid 
the  debugging  of  the  new  computer  code;  it  also  validated  the  new  computer  code  when 
the  two  results  matched  within  the  single-precision  NJOY  uses  for  calculations.  The 
interpolation  of  the  tabulated  data  and  other  intermediate  values  were  validated  using 
these  two  tools. 

IV.2.2:  NJOY  Validation  Using  Legendre  Moments 

The  second  validation  method  used  the  mechanism-specific  Legendre  moments 
output  that  is  generated  by  NJOY.  These  NJOY  Legendre  moments  could  then  be 
compared  to  two  different  computations  of  Legendre  moments  by  the  new  computer 
code.  The  new  computer  code  matched  the  Legendre  moments  from  NJOY  to  the 
practical  extent  possible. 

NJOY  has  several  behaviors  that  make  the  comparisons  difficult.  The  first 
NJOY  behavior  is  the  elimination  of  the  output  when  the  value  of  the  cross  section 
drops  below  10~9  barns.  This  behavior  creates  a  difficulty  because  some  mechanism  and 
group-to-group  pairs  have  values  below  the  NJOY  cutoff.  These  could  not  be  checked. 
Therefore,  I  was  forced  to  rely  upon  the  previous  validations  of  the  cross  section  curves 
and  interpolations,  and  the  validation  of  the  group-to-group  pairs  above  the  cutoff. 

Another  behavior  that  makes  comparison  to  NJOY  difficult  is  limited  error 
control  within  the  program.  For  example,  NJOY  will  perform  some  integrations  using 
a  pre-selected  (hard  coded)  number  of  intervals  for  a  composite  quadrature.  The  code 
does  not  refine  that  mesh  and  compare  quadrature  results  to  test  for  convergence,  nor 
does  it  use  any  other  error  control  scheme.  If  the  cross  section  involved  is  dominated 
by  other  mechanisms,  this  can  be  a  practical  approach,  but  it  does  not  provide  a  reliable 
benchmark.  In  such  cases,  the  NJOY  result  may  only  be  good  to  one  or  two  digits  (out 
of  the  four  digits  printed). 
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The  new  computer  code  was  designed  to  calculate  and  output  the  Legendre 
moments  for  comparison  to  the  NJOY  Legendre  moments.  The  comparisons  between 
NJOY  and  the  direct  calculations  for  the  Legendre  moments  from  my  code  were 
favorable  and  matched  to  within  0.1%.  I  considered  this  level  of  accuracy  to  validate  the 
integration  routines  and  the  joint  distributions  that  generate  the  Legendre  moments. 
Table  1  presents  a  summary  of  an  example  isotope,  mechanism,  and  energy  group-to- 
group  pair. 

Validating  the  integration  routines  that  generate  the  Legendre  moments  is  a 
necessary  condition  for  the  PAX  cross  sections  to  be  correct.  But,  I  did  not  consider  it 
sufficient.  The  approximation  to  the  Legendre  moments  using  PAX  cross  sections  from 
equation  (24)  presented  in  chapter  2,  were  also  generated  and  compared  to  the  NJOY 
moments. 


Method  of  Generation 

0th  moment 

1st  moment 

2nd  moment 

3rd  moment 

NJOY 

0.02432 

-0.01391 

0.002504 

0.001285 

New  code  using  direct 
calculation 

0.024317 

-0.013913 

0.0025041 

0.0012850 

New  code  with  PAXm 
approximation 

0.024317 

-0.013905 

0.0024859 

0.0013161 

Table  1:  Comparison  of  Legendre  moments 

10B  elastic  scatter  cross  sections  (in  barns)  for  group  1  to  4  using  LANL- 
30  defined  in  appendix  D 


All  of  the  ENDF/B-VI  representations  were  examined  with  results  similar  to 
this  example.  I  consider  the  new  computer  code,  including  the  new  PAX  cross  section 
approximation,  to  have  been  successfully  validated  using  the  combination  of  the 
different  comparisons. 

IV. 3:  Discrete  Elements  Validation 

Although  the  PAX  cross  sections  were  validated  by  the  investigation  using 
NJOY,  the  PAX  cross  section  used  to  generate  bin-to-bin  cross  sections  were  also 
validated.  I  chose  DelGrande’s  Monte  Carlo  cross  section  code  (7)  to  use  as  a 
benchmark  for  bin-to-bin  cross  sections  because  it  has  been  previously  validated  and 
published. 
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The  element-to-element  conditional  probability,  h  (tri),  was  also  validated 

directly  by  comparing  the  discrete  quadrature  one-dimensional  code  output  to  Monte 
Carlo  estimates  of  h  (rri)  for  all  (n  $n  )  pairs  of  a  slab-geometry  DEs  angular 

quadrature  partition  of  the  sphere  (14).  These  agreed  within  the  estimated  uncertainties 
of  the  Monte  Carlo  calculations. 

The  approximation  of  the  bin-to-bin  cross  sections  using  PAX  cross  sections 
presented  in  equation  (22)  agreed  within  the  estimated  uncertainties  of  DelGrande’s 
Monte  Carlo  cross  section  code.  I  consider  the  bin-to-bin  cross  sections  calculated 
using  the  PAX  cross  sections  to  be  valid  using  the  combination  of  the  results  from  the 
NJOY  comparison,  the  comparison  to  the  DelGrande’s  Monte  Carlo  cross  section  code, 
and  the  direct  comparison  of  the  element-to-element  conditional  probability  function  to 
Monte  Carlo  calculations. 
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V.  Performance  of  the  PAX  Cross  Section  Code 


The  performance  scaling  and  the  typical  memory  requirements  for  the  PAX 
cross  section  portion  of  the  scatter  cross  section  code  are  presented.  The  code  is 
examined  to  determine  how  the  computer  runtime  scales  with  respect  to:  the  number  of 
energy  groups,  the  number  of  tabulated  PAX  cross  section  points,  and  the  desired 
accuracy  of  the  calculation. 

The  runtime  scaling  power,  p,  can  be  empirically  estimated  from  two 
computational  runs  with  a  parameter  changing  from  N\  to  No  and  runtimes  D  t j  and 


Dt 


2 


hence, 


IVjf  = 
0 


Log(pt2/Dt j) 

Nl/Ni 


(44) 


(45) 


All  of  the  examples  in  the  following  sections,  with  the  exception  of  the 
investigated  parameter,  have  been  run  with:  the  30  group  structure  presented  in 
appendix  D,  a  material  temperature  of  300K,  a  relative  error  tolerance  of  0.00 1 ,  and  a 
PAX  mesh  of  64  equal-width  pieces.  The  computations  were  run  on  a  computer  with  a 
1  GHz  processor  and  512  MB  of  RAM.  The  examples  are  typical  of  other  testing  and 
evaluation  performed. 


V.  i :  Scaling  with  the  Number  of  Energy  Groups 

The  runtime  scales  as  approximately  linear  with  the  number  of  energy  groups. 
As  the  number  of  groups  increases  from  G,  to  Ga,  the  ratio  of  the  non-zero  group-to- 

group  pairs,  N^aus  to  N 2airs  ,  in  the  PAX  cross  section  should  vary  as 

pairs 
N  pairs 

This  does  not  result  in  quadratic  scaling  of  the  amount  of  computational  effort  because 
the  amount  of  computation  required  for  each  group-to-group  pair  should  be 


fJ 

iM 


(46) 
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proportional  to  ,  hence,  inversely  proportional  to  the  number  of  energy  groups,  G. 

The  offset  in  the  computational  effort  for  each  group-to-group  pair  combined  with  the 
increase  in  the  number  of  non-zero  group-to-group  pairs  results  in  the  observed 
approximately-linear  scaling 

As  an  example,  cross  sections  of  the  isotope  56Fe  were  generated  for  30,  59,  and 
117  groups  using  the  new  energy  groups  listed  in  appendix  D.  This  isotope  is 
representative  because  it  has  a  parameterized  resonance  region.  Additionally,  the  56Fe 
isotope  has  many  different  scatter  mechanisms.  Therefore,  the  runtime  scaling  is  a  good 
average  over  many  different  effects. 


Number  of  Groups 

Computational  Time 

Scaling  Power  (p) 

30 

3246  seconds 

59 

5193  seconds 

0.69 

117 

10079  seconds 

0.97 

Figure  8:  Scaling  of  code  runtime  with  the  number  of  energy  groups 

V.2:  Scaling  with  the  Number  of  PAX  Cross  Section  Pieces 

The  scaling  with  the  number  PAXk  cross  section  pieces,  or  scaling  with  K,  was 
approximately  quadratic.  An  increase  in  K  linearly  increases  the  computational  cost  of 
the  angular  integration  and  also  linearly  increases  the  computational  cost  of  the  incident 
energy  integration.  Nested,  the  two  linear  increases  are  quadratic  scaling. 

As  an  example,  cross  sections  for  were  generated  with  K  equal  to  64,  128, 

and  256.  This  isotope  is  a  good  example  of  runtime  scaling  as  K  increases  because 
elastic  scatter  is  the  only  mechanism  available  and  elastic  scatter  is  present  for  all 
materials  and  all  incident  energy  groups. 


Number  of  Points 

Computational  Time 

Scaling  Power  (p) 

64 

6  seconds 

128 

24  seconds 

2.0 

256 

117  seconds 

2.3 

Figure  9:  Scaling  of  the  code  runtime  with  PAX  cross  section  pieces 
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V.3:  Impact  of  Relative  Error  Tolerance  on  Runtime 


There  is  little  increase  in  the  runtime  with  lower  relative  tolerances.  I 
investigated  relative  tolerances  only  as  low  as  O.OOOl  because  the  ENDF/B-VI  data 
upon  which  the  PAX  cross  sections  depend  are  accurate  to  at  most  4  digits.  When  the 
integrations  are  performed  for  the  PAX  cross  sections,  the  relative  error  in  the 
integrations  usually  met  the  tightest,  0.0001,  relative  tolerance.  Therefore,  relaxing  the 
relative  tolerance  does  not  improve  runtime  drastically  because  often  no  extra  recursion 
was  required  to  pass  the  stricter  tolerance. 

As  an  example,  cross  sections  of  the  isotope  56Fe  were  generated  with  relative 
tolerance  settings  of  0.01,  0.001,  and  0.0001. 


Relative  Error  Tolerance 

Computational  Time 

0.01 

3014  seconds 

0.001 

3245  seconds 

0.0001 

4395  seconds 

Figure  10:  Runtime  versus  relative  error  tolerance 
V.4:  Memory  Requirements 

The  memory  requirements  for  the  PAX  cross  section  code  were  quite  modest. 
The  most  memory  required  for  any  of  the  materials  and  user-set  input  data  was  40  MB 
of  RAM.  Typically,  only  10  MB  were  needed.  For  practical  work,  the  memory 
requirements  are  not  limiting;  the  use  of  the  dynamically  allocated  arrays  presented  in 
chapter  3  was  successful. 
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VI.  1-D  Transport  Comparisons  for  Multi-Group  Cross  Section 

Approximations 

This  chapter  presents  two  test  problems  to  demonstrate  the  improvement  of 
using  the  discrete  elements  approximation  with  PAX  cross  sections  over  the  discrete 
ordinates  approximation  with  truncated  Legendre  expansions  for  the  group-to-group 
scatter  cross  sections.  One-dimensional  slab  geometry  transport  favors  the  discrete 
ordinates  approximation  in  comparison  to  the  discrete  elements  approximation.  The 
discrete  ordinates  Gauss-Legendre  angular  quadrature  should  approximate  the  angular 
flux  solution  better  than  the  discrete  elements  composite  midpoint  rule  because  the 
solution  for  y g  (x,  in)  is  continuous  in  x  and  m,  except  at  m  =  0 ,  given  incident  fluxes 
and  emission  sources  that  are  continuous  in  m.  Therefore,  if  discrete  elements  with 
PAX  cross  sections  produces  fluxes  that  are  at  least  as  good  as  those  produced  by 
discrete  ordinates  with  truncated  Legendre  group-to-group  scatter  cross  sections,  then 
I  anticipate  that  this  will  also  be  true  in  two-  and  three-dimensional  transport,  which 

favor  the  discrete  elements  approximation  because  yg  (r,W)  need  not  be  continuous  in 
r  and  W. 

The  two  test  problems  use  many  of  the  same  parameters  for  the  cross  sections 
and  the  transport. 

1.  Cross  sections  are  Doppler  broadened  to  300  K. 

2.  0.001  relative  tolerance  for  piecewise  average,  group  cross  sections 

3.  1.0'  10"  5  relative  tolerance  for  transport  convergence 

4.  P AX(54  cross  sections 

5.  Symmetry  boundary  at  left  end 

6.  Vacuum  boundary  at  right  end 

7.  Energy-weight  function,  F  (E)= 

8.  Isotropic  source  in  energy  group  1  of  the  new  30,  59,  and  117  group  structures 
defined  in  appendix  D 
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The  source  emits  only  in  group  1  so  that  the  computational  artifacts  of  the 
anisotropic  down-scatter  into  the  lower  energy  groups  would  not  be  obscured  by  a 
source  in  those  groups. 


VI.  i:  Test  Problem  i:  Thin  Source  Embedded  in  Water 


The  first  test  problem  has  water  throughout  the  entire  problem  with  a 
symmetry  boundary  on  the  left  and  a  vacuum  boundary  on  the  right.  The  source  is 
located  in  a  thin  region  of  water  on  the  left  side  of  the  problem  emitting  isotropically  in 
energy  group  1.  The  dimensions  in  the  figure  below  are  in  centimeters. 

Source 


Symmetry 


Vacuum 


0.0  0.1 


20.1 


Figure  1 1:  Thin  source  embedded  in  water 

Three  energy  groups  out  of  the  30-group  structure  were  chosen  for  comparison 
in  this  problem.  The  groups  were  chosen  to  demonstrate  the  differences  between  the 
anisotropic  and  isotropic  computations  as  well  as  the  differences  between  the  discrete 
ordinates  and  discrete  elements  computations.  Energy  group  1  was  chosen  to  show  the 
impact  of  the  anisotropy  when  the  scalar  flux  solution  is  dominated  by  the  removal 
cross  section  and  streaming.  Energy  group  2  was  chosen  to  show  the  impact  of  the 
anisotropy  of  the  down-scatter  source.  Finally,  energy  group  28  was  chosen  to  show 
that,  although  the  anisotropy  of  the  higher  energy  groups  is  important  to  the  value  of 
the  cross  section,  the  overall  down-scatter  source  is  mostly  isotropic  and  both  discrete 
ordinates  and  discrete  elements  perform  well. 

VI.  1.1:  Isotropic  Comparison 

The  first  series  of  figures  compares  calculations  using  both  isotropic  discrete 
ordinates  (S12/P0)  and  isotropic  discrete  elements  (DE12). 


44 
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Figure  14:  Isotropic  comparison  for  group  28 

The  good  agreement  between  both  the  discrete  elements  and  the  discrete 
ordinates  for  all  three  energy  groups  demonstrates  DEi2  is  an  adequate  resolution 
compared  to  the  discrete  ordinates  S12. 

VI.  1.2:  Anisotropic  Comparison 

The  anisotropic  comparison  includes  a  Monte  Carlo  multi-group  transport 
calculation  using  PAX64  cross  sections.  Energy  groups  1  and  2  exhibit  a  computational 
artifact  for  both  discrete  ordinates  and  discrete  elements  compared  to  the  Monte  Carlo 
solution  for  the  scalar  flux  due  to  the  optically  thin  source  region  (in  these  groups).  In 
one-dimensional  transport  with  a  thin  source,  much  of  the  contribution  to  the  scalar 
flux  in  and  near  the  source  region  arises  from  neutrons  traveling  nearly  perpendicular 
to  the  axis  ( |m|  =  1 ).  The  Monte  Carlo  simulation  includes  particles  emitted  with  \m\ 

near  or  equal  to  zero,  but  the  discrete  ordinates  (Sn  quadratures)  and  discrete  elements 
methods  need  very  high  angular  resolution  (many  elements  or  ordinates)  to  include 
ordinates  close  enough  to  perpendicular  near-perpendicular  contribution  to  the  scalar 
flux. 
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Instead  of  using  the  typical  discrete  ordinates,  Sn,  quadrature  based  on  a  Gauss- 
Legendre  quadrature  rule  (in  one-dimension),  a  two  angular  region  Gauss-Legendre 
rule  could  be  used.  The  discontinuity  in  the  angular  flux  at  m  =  0  can  be  better 
approximated  using  a  double  S„  quadrature  rule  (DS„),  such  as  the  DS6  quadrature  rule 
given  in  appendix  E.  The  DS8  quadrature  with  the  Ps  Legendre  expansion  for  the 
group  scatter  cross  section  has  been  included  to  examine  the  effects  of  a  quadrature  rule 
that  better  approximates  the  discontinuity  at  m  =  0  . 


x  (cm) 


Figure  15:  Energy  group  1  anisotropic  comparison 

The  first  energy  group  exhibits  the  angular  quadrature  artifact  in  the  scalar  flux, 
although  the  DS6  is  significantly  better  than  the  DEi2  or  the  any  of  the  Si2  quadratures. 
Near  the  source  region,  the  discrete  elements  and  the  discrete  ordinates  underestimate 
the  scalar  flux.  Note  that  the  DEJ2  is  a  better  approximation  than  the  Si2  quadratures. 

Although  the  DS6  quadrature  is  qualitatively  a  better  approximation  to  the 
scalar  flux  in  energy  group  1  than  DEi2,  an  angular  flux  comparison  demonstrates  the 
negativity  associated  with  the  P5  Legendre  expansion  for  the  group  scatter  cross 
section.  Figure  16  through  Figure  18  show  the  angular  flux  as  a  function  of  both  space 
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and  the  ordinate  or  element.  Note  that  the  angular  flux  is  an  angle  integrated  value  for 
the  discrete  elements  and  a  point  value  for  the  discrete  ordinates. 


Figure  16:  Energy  group  1,  DSs/Pb  positive  angular  flux 

The  white  areas  in  the  plot  are  negative  values  in  the  angular  flux  arising  from 
the  P 5  Legendre  expansion  for  the  group  scatter  cross  sections.  The  negative  values  for 
the  angular  flux  do  not  contribute  significantly  to  the  scalar  flux,  but  are  completely 
unphysical.  The  negative  angular  flux  is  plotted  in  Figure  17. 
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Figure  17:  Energy  group  1,  DSe/Pj  negative  angular  flux 

The  DE12  angular  quadrature  with  the  PAX64  cross  sections  result  in  positive 
angular  fluxes  throughout  the  entire  problem  when  a  source  is  present.  The  negative, 
unphysical  angular  flux  is  present  in  the  DS6/Ps  and  not  in  the  DEia. 
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Figure  18:  Energy  group  1,  DE12  angular  flux 

Energy  group  2  is  used  for  two  comparisons,  the  importance  of  the  anisotropy  to 
the  scalar  flux  solution  (thus  motivating  my  work  in  investigating  anisotropic  scatter 
cross  sections)  and  the  point-wise  convergence  of  the  scalar  flux  when  using  the 
truncated  Legendre  expansions.  The  down-scatter  source  from  group  1  is  highly 
anisotropic  and  the  scalar  flux  solution  is  incorrect  for  the  isotropic  approximation 
(shown  as  S12/P0). 
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Figure  19:  Energy  group  2  with  discrete  ordinates 

The  point-wise  convergence  of  the  scalar  flux  when  using  discrete  ordinates 
with  truncated  Legendre  expansions  is  evident  for  the  typical  Sn  quadrature.  The  DSe 
discrete  ordinate  quadrature  is  a  good  approximation  for  the  scalar  flux.  But,  the  S12 
quadratures  are  not  converging  uniformly  to  the  scalar  flux  solution  when  the 
truncation  order  of  the  Legendre  expansion  is  increasing.  The  P3  approximation  is 
closer  to  the  Monte  Carlo  scalar  flux  solution  for  spatial  positions  greater  than  1  cm  and 
the  P5  approximation  is  closer  for  spatial  positions  less  than  1  cm.  The  Po  and  Pi 
approximations  do  not  show  any  uniform  convergence  to  the  scalar  flux  compared  to 
the  other  truncated  Legendre  expansions.  Thus,  different  discrete  ordinates 
quadratures  (S12  and  DSe)  combined  with  different  truncated  Legendre  expansions  do 
not  converge  uniformly  to  the  scalar  flux  solutions — leading  to  the  art  form  of  choosing 
an  appropriate  truncation  order  and  discrete  ordinate  quadrature  set. 

Having  shown  the  non-uniform  convergence  to  the  Monte  Carlo  solution  for  the 
scalar  flux,  an  investigation  of  the  convergence  was  made  for  the  discrete  elements 
approximations  using  the  PAX64  cross  sections.  The  scalar  flux  for  five  discrete 
elements  quadratures  was  compared  to  the  scalar  flux  estimate  from  Monte  Carlo  in 
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energy  group  2.  The  discrete  elements  angular  quadratures  are  converging  towards  the 
Monte  Carlo  estimate  for  the  scalar  flux. 


Figure  20:  Energy  group  2  with  discrete  elements 

Energy  group  28  shows  the  continued  importance  of  the  anisotropic  cross 
sections.  Although  the  Pi  expansion  is  different  from  the  other  calculations,  the  discrete 
elements  and  discrete  ordinates  calculations  match  the  Monte  Carlo  estimate  of  the 
scalar  flux  closely. 
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Figure  21:  Energy  group  28  anisotropic  comparison 

The  thin  source  in  a  water-medium  test  problem  displayed  two  important 
features  of  the  different  discrete  approximations  to  the  scalar  flux  solution.  The 
anisotropic  scatter  cross  section  dominates  the  shape  of  the  scalar  flux.  And,  the 
discrete  elements  approximation  to  the  scalar  flux  is  converging  to  the  Monte  Carlo 
solution  while  the  discrete  ordinates  to  the  scalar  flux  using  truncated  Legendre 
expansions  is  not  converging  uniformly. 

VI. 2:  Test  Problem  2:  Multi-Layer  Shield 

Five  different  effects  were  examined  with  the  second  test  problem:  negative 
scalar  fluxes  when  using  diamond  difference  spatial  quadratures,  incorrect  scalar  fluxes 
when  using  discrete  ordinates  with  truncated  Legendre  expansions  for  the  scatter  cross 
sections,  convergence  to  the  wrong  energy  dependence  when  using  truncated  Legendre 
expansions,  uniform  convergence  of  the  energy  dependence  with  increased  number  of 
discrete  elements,  and  the  importance  of  low  probability  scatters  to  the  energy 
dependence. 
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This  problem  used  four  materials  as  a  shield  for  a  high  energy  source.  Figure  22 
shows  the  diagram  for  test  problem  2  where  the  dimensions  are  all  in  centimeters.  The 
isotropic  source  in  energy  group  1  was  distributed  uniformly  throughout  the  56Fe 
region.  All  of  the  examples  used  96  equal-thickness  spatial  cells. 

56Fe  ‘H2160  10B4C  207Pb 


Symmetry 


0.0  f  1.25 

Distributed 

Source 


_ Vacuum 

3.75  5.0  6.25 


Figure  22:  Multi-layer  shield 

VI. 2.1:  Diamond  Difference  Generates  Negative  Scalar  Fluxes 

The  first  effect  examined  was  the  negative  scalar  fluxes  when  the  diamond 
difference  (DD)  spatial  quadrature  method  was  used.  The  scalar  flux  for  an  exponential 
characteristic  (EC)  spatial  quadrature  method  is  compared  to  the  DD  calculation  for 
energy  group  30  with  both  spatial  quadrature  methods  using  a  discrete  element  angular 
quadrature  with  PAX61,  cross  sections  (ie.  non-negative  cross  sections). 
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Figure  23:  EC  comparison  to  DD 

The  DD  calculation  is  incorrect  throughout  the  entire  problem  for  this  energy 
group.  The  scalar  flux  in  the  B4C  region  (between  3.75  and  5.0  cm)  is  unphysical  and 
oscillates  between  positive  and  negative  values,  as  is  seen  in  Figure  24  (an  enlargement 
of  the  B4C  region). 

The  oscillations  of  the  DD  method  in  the  B4C  region  have  led  to  a  drastic 
overestimation  of  8  orders  of  magnitude  in  the  scalar  flux  in  the  207Pb  region  (between 
5.0  to  6.25  cm).  The  DD  method  has  failed  to  produce  positive  scalar  fluxes  even  with 
non-negative  group-to-group  cross  sections. 
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Figure  24:  DD  close-up  of  B4C 

The  test  problem  was  run  with  continued  spatial  refinements  using  DD  to 
discover  when  the  scalar  flux  in  energy  group  30  would  no  longer  be  negative.  The  DD 
scalar  flux  solution  continued  to  have  negative  values  for  the  scalar  flux  when  over  5000 
spatial  cells  were  used — a  physical  thickness  of  just  over  10  mm  per  cell.  Higher 
refinements  were  not  possible  with  the  spatial  transport  code  because  the  code  was  not 
written  for  efficient  memory  usage,  but  rather,  it  stores  all  variables  to  facilitate 
examination  of  them. 

VI. 2. 2:  Discrete  Ordinates  Compared  to  Discrete  Elements 

The  step  spatial  quadrature  method  was  used  to  compare  discrete  ordinates  with 
truncated  Legendre  expansions  for  the  cross  sections  to  discrete  elements  using  PAX64 
cross  sections.  The  step  method  was  used  because  it  is  a  non-negative  method  (as 
opposed  to  DD)  with  non-negative  cross  sections.  An  EC  calculation  was  also 
performed  for  comparison.  The  EC  method  is  also  non-negative,  but  the  EC  method 
cannot  use  negative  cross  sections  and  the  step  method  can. 
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The  multi-layer  shield  problem  was  run  using  the  step  spatial  quadrature  for 
both  S12/P5  and  DE  12.  Energy  group  30  was  again  examined  and  the  results  are  shown 
in  Figure  25. 


Figure  25:  Step  comparison  to  EC 

The  scalar  flux  for  the  step  method  is  close  to  the  EC  scalar  flux  solution  when 
both  methods  used  the  discrete  elements  angular  quadrature.  This  indicates  the 
inadequacy  of  the  spatial  mesh  for  the  0(s  Dr)  step  method  as  compared  with  the 

0(s  Dj)4  EC  method.  By  comparison,  the  much  larger  difference  between  the 

step/DE  12  and  the  step/Sia/Ps  results  is  significant.  I  attribute  this  difference  to  the 
inaccuracies  of  the  discrete  ordinates  angular  quadrature  with  a  Ps  truncated  Legendre 
expansion. 
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Figure  26:  Close-up  of  step  comparison  to  EC 

An  enlargement  of  the  scalar  flux  solution  for  the  B  tC  and  207Pb  regions  of  the 
scalar  flux  shows  the  inaccurate  scalar  flux  using  the  S12  angular  quadrature  with  P 5 
truncated  Legendre  expansion  in  the  outer  regions. 


VI. 2. 3:  Effects  of  Angular  Quadrature  Refinement  on  Energy  Dependence 

The  effects  of  the  angular  quadrature  refinement  on  the  energy  dependence  were 
examined  by  comparing  the  energy  dependence  using  a  step/DEia  to  four  step/Sn/Pi 
calculations.  The  EC/DE12  agreed  with  the  step/DEia  calculation.  The  step  spatial 
quadrature  method  is  again  used  because  of  its  non-negativity  when  using  non-negative 
cross  sections.  The  30-group  structure  from  appendix  D  is  used  in  Figure  27  to  plot  the 
neutron  current  exiting  the  test  problem  at  the  vacuum  boundary. 
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Figure  27:  Effect  of  discrete  ordinates  angular  resolution  in  energy  dependence 

The  angular  quadrature  is  not  converging  to  the  step/DEia  solution  for  the  four 
step/Sn/Pi  calculations.  This  is  shown  by  the  more  accurate  S2/P1  solution  for  energies 
less  than  10  eV  and  the  less  accurate  S2/P1  solution  for  energies  greater  than  10  eV.  To 
investigate  if  the  discrete  ordinates  solutions  with  the  truncated  Legendre  expansions 
were  consistently  differing  from  the  discrete  elements  solutions,  a  plot  of  the  relative 
error  e  ^  of  the  discrete  ordinates  compared  to  the  discrete  elements  was  made  in 
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Figure  28:  Relative  difference  between  step/Sn/Pi  and  step/DEia 

The  relative  difference  between  the  step  with  different  Sn/Pi  combinations  and 
the  step  spatial  quadrature  with  DE12  shows  that  some  combinations  are  less  accurate  at 
some  energies  and  more  accurate  at  other  energies.  The  inclusion  of  the  double  Sn 
quadrature  (DS12)  shows  that  the  artifact  discussed  in  test  problem  1  is  not  dominating 
the  error  in  the  solution.  Instead  the  difference  between  the  step/DEis  and  the 
step/Sn/Pi  can  be  attributed  to  the  truncated  Legendre  expansions  for  the  group  scatter 
cross  sections. 


VI. 2.4:  Angular  Refinement  of  Discrete  Elements  with  Energy  Dependence 

Convergence  of  the  discrete  elements  angular  quadrature  is  expected  as  more 
elements  are  added.  For  this  test  problem,  Figure  29  shows  the  EC  spatial  quadrature 
methods  with  DE2,  DE4,  DE@,  and  DE12. 
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Figure  29:  Effect  of  discrete  elements  angular  resolution  in  energy  dependence 

In  a  similar  comparison  to  Figure  28,  relative  differences  between  the  DEa,  DE4, 
and  DE@  to  the  DE  12  solution  are  shown  in  Figure  30. 
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Figure  30:  Relative  difference  between  DEi2  and  DE2,  DE4,  and  DE6 

The  error  in  the  DE2  calculation  is  dominated  by  its  nearly  isotropic 
approximation.  The  angular  quadrature  for  DE2  is  insufficiently  angularly  resolved  for 
this  test  problem.  The  DE4  and  the  DE6  are  converging  to  the  DEi2  solution  for  the 
energy  dependence  of  the  current  at  the  vacuum  boundary. 


VI. 2. 5:  Energy  Group  Refinement 

The  effects  of  refining  the  energy  groups  for  this  test  problem  were  investigated. 
Three  different  comparisons  with  energy  group  refinement  were  made:  EC  with  DEi2, 
step  with  S12/P5,  and  step  with  DEi2.  The  energy  meshes,  with  30,  59,  and  117  groups, 
were  made  using  the  energy  group  structures  defined  in  appendix  D. 
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Figure  33:  Energy  dependence  for  step  with  DEia 

All  three  energy  dependencies  are  changing  with  the  refinement  in  energy 
group.  Although  the  step/Sia/Ps  appears  qualitatively  similar,  the  comparison  of 
step/Si2/Ps  to  the  step/DEia  in  section  VI. 2. 3  demonstrated  that  the  current  was 
incorrect.  The  energy  dependence  is  changing  for  both  the  step/DEia  and  EC/DEia 
calculations,  as  the  energy  group  structure  is  refined.  Therefore,  the  energy  group 
structure  may  not  be  sufficiently  refined  in  energy  to  accurately  determine  the  energy 
dependence  of  the  current  at  the  vacuum  boundary,  depending  on  the  application. 


VI. 2. 6:  Attributing  Features  in  the  Energy  Dependence 

A  comparison  was  made  between  the  rightward  partial  current  entering  the 
2°7pb  region  and  the  rightward  partial  current  at  the  vacuum  boundary  because  the 
refinement  in  energy  in  Figure  31  showed  an  increasingly  variable  solution  in  the 
energy  dependence  in  the  energy  range  of  10  keV  to  1  MeV.  The  comparison  was  made 
using  the  1 17-group  structure  presented  in  appendix  D.  Figure  34  shows  the  rightward 
current  at  both  the  vacuum  boundary  and  the  Bi,C/207Pb  boundary. 
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Figure  34:  Comparison  of  energy  dependence  of  rightward  partial  current  at  two 
spatial  locations 

Without  any  intrinsic  source  in  the  207Pb  region,  the  variation  of  the  energy 
dependence  of  the  current  exiting  at  the  vacuum  boundary  arises  from  either  the  cross 
section  of  207Pb  or  the  current  flow  from  the  B4C.  Figure  34  shows  that  the  current 
entering  the  207Pb  from  the  B4C  is  smooth  in  the  10  keV  to  1  MeV  energy  range  in 
comparison  to  the  current  at  the  vacuum  boundary.  Due  to  the  smoothness  of  the 
entering  flux  and  the  large  variation  of  the  total  cross  section  of  207Pb  shown  in  Figure 
35  in  this  energy  range,  I  attribute  the  variation  in  the  energy  dependence  of  the 
rightward  partial  current  at  the  vacuum  boundary  to  the  207Pb  cross  section. 
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Figure  35:  207Pb  total  cross  section 

Figure  34  shows  that  energy  dependence  for  energies  less  than  100  eV  of  the 
rightward  partial  current  at  the  vacuum  boundary  is  fairly  flat  compared  to  the  current 
exiting  the  B  iC  region.  Because  the  current,  in  this  problem,  exiting  at  the  vacuum 
boundary  can  only  arise  from  inflow  from  the  BiC  region  or  down-scatter  within  the 
207Pb  region,  the  low  energy  current  must  be  down-scatter  because  the  inflow  current  is 
negligible. 

To  examine  the  contributions  to  the  down-scatter  source,  a  calculation  with  the 
30-group  structure  in  appendix  D.2  was  used.  The  30-group  calculation  has  similar,  but 
less  resolved  energy  dependence  of  the  currents.  The  207Pb  region  down-scatter 
contribution  to  group  30  from  each  incident  group  was  calculated  as  an  angle- 
integrated,  region-integrated  down-scatter  source.  To  demonstrate  the  low  probability 
of  scatter  from  the  incident  group,  the  fractions  of  the  scatter  cross  section  from  the 
incident  group  to  group  30  are  also  presented  in  Table  2. 

Table  2:  Down-scatter  source  contribution  to  207Pb  region 
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Incident  Group 

Down-Scatter  Source  Fraction  of  Scatter  from 

Incident  Group  to  Group  30 

1 

6.755617E-08 

1.274419E-09 

2 

2.545  lllE-08 

1.320084E-09 

3 

1.559567E-08 

1.40671  lE-09 

4 

1.341229E-08 

2.325594E-09 

5 

9.901283E-09 

2.036932E-09 

7 

9.066829E-14 

9.475451E-15 

8 

6.584222E-13 

1.225620E-13 

9 

1.137028E-13 

1.919159E-14 

29 

4.874568E-10 

9.697575E-03 

The  group-29  down-scatter  source  was  smaller  than  the  down-scatter  source  for 
each  incident  group  1  through  5.  Therefore,  I  attribute  the  energy  dependence  of  the 
current  below  100  eV  to  the  down-scatter  source  from  energies  in  the  range  of  7.79  to 
14.1  MeV.  Low-probability  scatters  at  the  high  energies  are  dominating  the  energy 
dependence  at  low  energies  of  the  current  at  the  vacuum  boundary.  Among  109  group- 
one  neutrons  that  scatter  in  207Pb,  only  1.3  neutrons  scatter  directly  into  group  30. 
Nevertheless,  in  this  problem,  more  than  half  (51%)  of  the  neutrons  that  enter  group  30 
in  the  lead  do  so  by  direct  scatter  from  group  1. 
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VII.  Summary  and  Conclusions 


Nonlinear,  accurate,  and  robust  spatial  quadrature  schemes  for  the  multi-group 
discrete  ordinates  method,  such  as  the  exponential  characteristic  method,  provide 
physically  meaningful,  non-negative  fluxes  given  non-negative  incident  fluxes,  emission 
sources,  and  scattering  cross  sections.  Conventional  techniques  based  on  Legendre 
expansions  of  the  scattering  cross  sections  cannot  be  relied  upon  to  provide  such  cross 
sections.  Multi-group  discrete  element  bin-to-bin  scattering  cross  sections,  introduced 
by  DelGrande  and  Mathews  (7),  are  non-negative  if  calculated  from  non-negative 
scattering  cross  sectionss5  ( E E,m ).  These  cross  sections  eliminate,  or  at  least 

strongly  ameliorate,  many  other  artifacts  of  conventional  discrete  ordinates  cross 
sections  and  angular  quadratures,  as  discussed  by  DelGrande  and  Mathews  (7). 

The  original  objective  of  this  effort  was  to  develop  an  accurate  and  efficient 
scheme  for  calculating  the  cross  sections  directly  from  the  data  in  ENDF/B-VI  (18). 
This  would  make  it  possible  to  use  the  exponential  characteristic  method,  and  similar 
methods,  for  anisotropic  scattering  in  real  materials.  This  has  been  accomplished. 

DelGrande  and  Mathews  used  Monte  Carlo  simulations  that  sampled  from  ACE 
files  (13).  I  explored  enhancing  his  approach  by  introducing  variance  reduction 
schemes,  but  abandoned  that  approach.  It  is  easy  to  obtain  modest  accuracy  in  the  cross 
sections  for  the  likely  scatters,  but  to  do  so  for  one-in-a-million  scatters  seems 
unachievable.  The  multi-layer  shield  problem  presented  here  shows  that  such  unlikely 
scattering  events  can  be  the  dominant  source  of  particles  in  at  least  some  energy  groups 
in  some  problems,  so  I  sought  a  method  that  would  calculate  such  cross  sections 
accurately  and  efficiently. 

By  separating  the  angular  and  energy  dependences  of  the  discrete  elements  and 
energy  group  structure,  I  found  it  possible  to  formulate  the  cross  sections  such  that  the 
maximum  nesting  of  integrations  was  only  three  deep,  and  usually  only  two.  Nested 
adaptive  numerical  quadratures  have  proven  effective  in  making  these  calculations 
efficient  while  providing  user-set  convergence  tolerances.  This  has  been  implemented 
in  FORTRAN-95,  demonstrating  the  computational  practicality  of  the  approach.  An 
efficient  code  for  generating  the  scattering-cosine-dependent  conditional  element-to- 
element  scattering  probabilities,  h  (ni),  is  the  only  piece  needed  for  efficient 
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generation  of  bin-to-bin  cross  sections  for  multi-dimensional  transport  calculations. 
Many  approaches  to  defining  partitions  of  the  unit  sphere  into  discrete  elements  are 
possible;  each  would  need  its  own  code  for  the  element-to-element  conditional 
probabilities. 

In  order  to  validate  my  computational  techniques  and  the  code  that  implemented 
them,  it  was  necessary  to  generate  Legendre  coefficients  for  group-to-group  cross 
sections  and  compare  them  with  the  values  produced  by  NJOY. 

By  introducing  a  scattering  operator,  it  became  clear  that  one  code  could  provide 
not  only  bin-to-bin  cross  sections  and  Legendre  moments,  but  also  group-to-group 
cross  sections  approximated  as  piecewise-average  functions  of  the  scattering  cosine, 
which  I  call  PAX  cross  sections.  The  PAX  cross  sections  are  one  of  the  major 
contributions  of  this  work.  They  are  readily  computed  and  tabulated  and  can  then  be 
used  to  produce  bin-to-bin  cross  sections  and  Legendre  moments.  Furthermore,  the 
PAX  cross  sections  could  be  used  directly  for  sampling  scattering  in  multi-group  Monte 
Carlo  transport  simulations,  or  could  be  used  to  generate  equal-likelihood  scattering 
intervals  for  use  in  such  codes. 

The  code  written  for  this  research  is  not  a  production  code.  Some  of  the 
ENDF/B-VI  scatter  laws  are  not  yet  implemented.  Also,  writing  independent  modules 
for  computing  Legendre  moments,  bin-to-bin  cross  sections,  and  piecewise-average 
cross  sections  could  substantially  reduce  run  times.  Nevertheless,  the  practicality  and 
some  of  the  benefits  of  my  methods  have  been  demonstrated  successfully  using  this 
code. 

The  impact  of  these  new  approaches  to  generating  and  using  scattering  cross 
sections,  if  adopted  by  the  community,  will  be  to  provide  physically  meaningful  and 
more  accurate  cross  sections  and  physically  meaningful  and  more  accurate  fluxes  from 
radiation  transport  calculations  that  use  them. 
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Appendix  A:  Variable  Definitions 

A.  1 :  Continuous  Representations 

y(r,E,ty  (48) 

The  angular  flux  is  a  distribution  function  in  space,  energy,  and  angle.  It  can  be 

neutron  cm 

described  as  a  neutron  path  length  rate  density  with  units  such  as  - ~ - . 

s cm  i  eV  steradian 

ss  (r,E£®  E,V$xW)  (49) 

The  continuous  scatter  cross  section  is  a  point  function  in  space,  incident  energy,  and 
incident  direction.  It  is  distribution  function  with  respect  to  the  secondary  energy  and 
the  secondary  direction.  The  continuous  scatter  cross  section  has  units  of 
barns 

eV  steradian 


A. 2:  Energy  Group  Representations 

y  (r,E,  W)»  yg(r, W)F(£),  Eg£  E£  Eg_x  (50) 

where  the  energy-dependent  spectral  weighting  function  is  normalized  by 


5  dE  F  (E)=  1 


(51) 


D  E„ 


and 


-S  dE  —  q  dE 


O 


(52) 


D  E„  E  . 

g  g-  l 


The  group  angular  flux,  y  „  ,  is  a  distribution  function  in  space  and  angle,  and  a  bin 


integrated  function  in  energy.  It  has  units  of 


neutron  cm 


s  cm 
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(53) 


5  dE(>  5  dEss  (r,E  0®  E,  \\$><W)F  (£  0) 


(?,V$xW)  = 


DE 


DE0 


5  dEfr(Efy 


D  E 


gt* 


The  group  scatter  cross  section  is  a  point  function  in  space  and  incident  direction,  and  a 
spectrum-weighted  average  over  energy.  It  is  a  distribution  function  with  respect  to  the 
secondary  direction  and  a  bin  integrated  function  with  respect  to  the  secondary  energy. 

barns 


The  group  scatter  cross  section  lias  units  of 


steradian 


A. 3:  Discrete  Element  Representations 


ynjh 

d  w„ 


(54) 


where  DW[;  is  an  element  of  solid  angle.  The  discrete  element,  group  angular  flux  yn  g 

is  a  distribution  function  in  space  and  a  bin-integrated  function  in  energy  and  angle.  It 
neutron  cm 

has  dimensions  of - ~ - . 

s  cm  3 


5  dW  q  dV#  q  dEtfr  dE  ss  (r,E  0®  E,V#X W)F(£0) 


r  dwj, 

s'  ,  (/')=  — - 

n%n,g%g  v  ' 


DW 


DE 


gt 


DE„ 


DW.6  dEfE(Et) 


(55) 


DE 


gt 


The  bin-to-bin  cross  section  is  a  point  function  in  space,  a  spectrum-weighted  average 
in  energy,  and  an  element  average  in  incident  direction.  It  is  a  bin-integrated  function 
with  respect  to  the  secondary  energy  and  the  secondary  direction.  The  bin-to-bin  cross 
section  has  units  of  barns. 
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Appendix  B:  Summary  of  FORTRAN  90/95  Syntax  Terms 

A  more  complete  presentation  of  FORTRAN  90/ 95  syntax  can  be  found  in  reference 

(7). 

•  Al  1  oca  table  Array 

An  array  whose  shape  and  size  are  not  determined  until  space  is  created  for  the 
array  by  means  of  an  allocate 

•  Alloca te(  array( ) ) 

Command  the  creates  an  array  at  execution 

•  Call  mysub( ) 

A  call  to  user-defined  subroutine,  rnysub 

•  Close(  my  file  ) 

A  command  to  close  a  file  previously  opened  for  1/ O 

•  Cycl e 

A  command  to  execute  the  next  loop  in  a  Do  constuct 

•  Do...End  Do 

A  loop  initiated  by  Do  and  terminated  by  End  Do 

•  Derived  Type 

A  programmer-defined  variable  type,  typically  a  compound  of  other  intrinsic 
types 

•  Extent 

Number  of  elements  in  a  particular  dimension  of  an  array 

•  If  (  )  Then...ElseIf  (  )  ...Else...End  If 
A  logical  construct 

•  Open(  myfile  ) 

A  command  to  open  a  file  for  I/O 

•  Read( ) 

A  command  to  retrieve  information  from  a  file  previously  opened  for  input 

•  Select  Case  (  )  ...Case  (  )  ...End  Select  Case  (  ) 

A  logical  construct  to  select  from  many  cases  using  the  same  variable 
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Appendix  C:  Doppler  Broadening 


The  cross  section  equations  presented  in  chapter  2  are  general  enough  to  include 

Doppler  broadening.  The  cross  section,  s  (E  0),  in  equation  (30)  changes  to  s  ( E  tfT  ), 

an  effective  cross  section  that  includes  the  temperature  dependence.  The  broadening  is 
important  to  the  cross  section  calculation  because  the  total  energy  available  for  the 
reaction  is  the  important  parameter  for  determining  the  likelihood  of  a  reaction. 
Therefore,  the  energy  of  the  incident  neutron  and  the  average  energy  of  the  target 
nuclei,  characterized  by  the  material  temperature,  combine  to  form  a  broadened  cross 
section  or  an  effective  cross  section.  This  broadening  of  the  cross  section  is  typically 
called  Doppler  broadening.  If  a  transport  calculation  is  to  be  performed  at  realistic 
temperatures  such  as  between  three  hundred  to  a  few  thousand  degrees  Kelvin,  then  the 
cross  section  must  be  broadened  in  order  for  the  transport  result  to  be  meaningful. 
Additionally,  this  research  required  validation  from  sources  that  contained  only  Doppler 
broadened  data. 

Prior  to  deriving  the  broadening  integral,  I  assume  that  the  functional 
dependence  of  the  secondary  neutron  energy  and  the  angle  of  scatter  of  the  neutron  do 
not  strongly  depend  on  the  temperature  of  the  interacting  material.  This 
approximation  exactly  follows  the  NJOY  development  (13).  The  value  of  the  cross 
section  is  the  only  function  broadened  and  the  angular  and  secondary  energy 
distributions  remain  unaffected  by  the  Doppler  broadening. 

The  effective  cross  section  for  a  material  at  temperature  T  is  defined  to  be  that 
cross  section  that  gives  the  same  reaction  rate  for  a  stationary  target  nuclei  as  the  real 
cross  section  gives  for  moving  nuclei.  The  effective  cross  section  is 

rv's (v,T  )=  q  dvir  |v  -  vj\s  (]v  -  v>  ( v<kT  )  (56) 

where  v  is  the  velocity  of  the  incident  neutron,  v  t  is  the  velocity  of  the  target,  r  is  the 

j- 

density  of  the  target  nuclei,  s  is  the  cross  section  for  stationary  nuclei,  and  P  (\>%T  )  is 

the  distribution  function  for  the  target  nuclei  velocities  in  the  laboratory  system. 

The  velocity  distribution  of  the  target  nuclei  can  be  a  very  complicated  function. 
For  many  applications,  however,  the  target  motion  can  be  approximated  as  isotropic 
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and  the  distribution  of  the  velocities  can  be  approximated  by  the  Maxwell-Boltzmann 
function, 

r  r  „3/  »  r 

P(v$T)dvt=  — av*dv$  (57) 

p 

M 

where  a  =  - ,  k  is  Boltzmann’s  constant,  and  Mbs  the  target  mass. 

2 kT  8 

Equation  (56)  can  be  partially  integrated  in  terms  of  the  relative  speed  defined  as 

V  =  |v 0-  v|,  (58) 

to  give  the  standard  Doppler  broadened  cross  section  as  the  sum 

■S' ( y,r  )=  5  *  (v,T  )-  s*  (-  v,T  ),  (59) 

with  the  function  defined  as 


*>,0=  ^4  dvstryh-^-^. 

p/' 2v  0 


(60) 


The  exponential  function  in  equation  (60)  limits  the  significant  part  of  the  integral  to 
the  range 


4  4 

v  -  - <  V  <  v  +  - . 

Vo"  s/a~ 


For  s  C-  v>,  the  integral  depends  only  on  velocities  satisfying 

4 


(61) 


0£  V  < 


■sfa 


(62) 


Because  all  of  the  data  used  in  the  ENDF/B-VI  are  given  energy,  rather  than  velocity, 
these  results  can  be  converted  to  energy  units  using 


1  s4o2  16 kT 

—m 


A 


(63) 


The  numerical  evaluation  of  equation  (60)  developed  for  NJOY  and  used  in  this 
dissertation  assumes  that  the  cross  section  can  be  represented  by  a  piecewise  linear 
function  of  energy  to  acceptable  accuracy.  A  piecewise  linear  function  to  represent  the 
resonance  regions  was  generated  to  satisfy  this  need.  The  method  used  to  generate  the 
piecewise  linear  function  representation  was  a  simple  bisection  method.  Having 
generated  a  piecewise  linear  representation,  equation  (60)  was  evaluated  using  the 
following  procedure. 
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Using  the  subscripts  in  the  following  derivations  to  denote  the  piecewise  linear 
values  obtained  for  the  resonance  regions  and  defining  the  reduced  variables  y  =  -Jav 
and  x  =  Vo" v  ,  the  cross  section  becomes 


with  slope 


5  00  =  s.  +  Si(x2-  xj), 
s  =  (!W  -  jl) 


Equation  (60)  can  now  be  written  as 


N  Xi+ 1 


1  iv  7 

^(v)=  — ra  5  dxs  ODxV  ix' yJ  =  a  {*.  | 
'Jpy  i=  o  v  i 


.  -  S.X. 

l  l  l 


,2}j+  Bf  },  (66) 


where 


x0  =  0  («■ 

+  1  =  ¥  •  (6I 

4  =  +  (6! 

«,  =  U/4  +  iff,  +  6tf ,  +  4v//,  +  y2H„,  (7< 

and  where  Hn  is  shorthand  for  Hn  (x.  -  -  _y\  The  extrapolations  to  zero  and 

infinity  assume  a  constant  cross  section  (sQ  =  s  ^  =  0).  The  H  functions  are  the 
incomplete  probability  integrals  defined  by 

Hn(g,b')=  ~^—r\dzzne'z  .  (7: 

Vp 

These  functions  can  be  computed  by 


where 


Hn  (a,b)=  F„  da|)-  Fn  (|&|> 


1  2 

Fn(ay=  - r\dzznez  . 

Jp 


The  Fn  functions  satisfy  a  recursion  relation  that  can  be  used  to  obtain 
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Fn  to >  =  —  erfcto), 

0  2 

F,  to>=  — - — exp  (-  a2\ 

1  Up  FV  J 

F„  to)  =  - - -F  0  to>+  a"'  1F1  to), 

n  ^  n  -  2  1  7 


and  where  the  erfc  to  )  denotes  the  complementary  error  function 

2  ¥  72 

erfc  to  )  =  - c\dze  "  . 


(74) 

(75) 

(76) 


(77) 
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Appendix  D:  Energy  Group  Structures 

All  of  the  energy  group  structures  used  throughout  the  dissertation  are 
documented  here.  The  values  are  all  in  eV  and  are  displayed  using  the  FORTRAN 
exponential  notation.  The  values  displayed  are  the  energy  group  boundaries. 


D.l:  LANL-30 


1.70E+7 

3.68E+6 

3.03E+5 

4.54E+2 

4.14E-1 

1.50E+7 

2.87E+6 

1.84E+5 

1.67E+2 

1.52E-1 

1.35E+7 

2.23E+6 

6.76E+4 

6.14E+1 

1.39E-4 

1.20E+7 

1.74E+6 

2.48E+4 

2.26E+1 

1.00E+7 

1.35E+6 

9.12E+3 

8.32EOO 

7.79E+6 

8.23E+5 

3.35E+3 

3.06E00 

6.07E+6 

5.00E+5 

1.24E+3 

1.13E00 

D.2:  Altered  LANL-30 

This  energy  group  structure  is  identical  to  the  LANL-30  shown  above  except 
that  the  highest  energy  groups  has  been  replaced  with  a  group  containing  14.0  to  14.1 
MeV  neutrons.  The  second  energy  group  has  also  changed  to  accommodate  the  new 
highest  energy  group. 


1.41E+7 

3.68E+6 

3.03E+5 

4.54E+2 

4.14E-1 

1.40E+7 

2.87E+6 

1.84E+5 

1.67E+2 

1.52E-1 

1.35E+7 

2.23E+6 

6.76E+4 

6.14E+1 

1.39E-4 

1.20E+7 

1.74E+6 

2.48E+4 

2.26E+1 

1.00E+7 

1.35E+6 

9.12E+3 

8.32EOO 

7.79E+6 

8.23E+5 

3.35E+3 

3.06E00 

6.07E+6 

5.00E+5 

1.24E+3 

1.13E00 

D.3:  New  59  Group  Structure 

This  group  structure  is  based  on  the  altered  LANL-30  presented  above  except 
that  each  energy  group  has  been  cut  into  two  pieces  using  the  geometric  mean.  The 
formula  is  the  square  root  of  the  product  of  the  energy  group  boundaries.  The  highest 
energy  group  has  not  been  cut  into  two  pieces. 
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1.41E+7 

3.68E+6 

3.03E+5 

4.54E+2 

4.14E-1 

3.249862E+6 

2.361 186E+5 

2.753507E+2 

2.508545E-1 

1.40E+7 

2.87E+6 

1.84E+5 

1.67E+2 

1.52E-1 

1.S7477SE+7 

2.529842E+6 

1.115276E+5 

1.012610E+2 

4.596520E-3 

1.S5E+7 

2.23E+6 

6.76E+4 

6.14E+1 

1.39E-4 

1.272792E+7 

1.969822E+6 

4.094484E+4 

3.725104E+1 

1.20E+7 

1.74E+6 

2.48E+4 

2.26E+1 

1.095445E+7 

1.532645E+6 

1.503915E+4 

1.371248E+1 

l.OOE+7 

1.S5E+6 

9.12E+3 

8.32EOO 

8.826098E+6 

1.054064E+6 

5.527386E+3 

5.04571  lE+OO 

7.79E+6 

8.23E+5 

3.35E+3 

3.06E00 

6.876431E+6 

6.414827E+5 

2.038 136E+3 

1.859516E+00 

6.07E+6 

5.00E+5 

1.24E+3 

1.13E00 

4.726267E+6 

3.892300E+5 

7.503066E+2 

6.839737E-1 

D.4:  New  i  17  Group  Structure 

This  group  structure  is  based  on  the  new  59  group  structure  presented  above 
except  that  each  energy  group  lias  been  cut  into  two  pieces  using  the  geometric  mean. 
The  formula  is  the  square  root  of  the  product  of  the  energy  group  boundaries.  The 
highest  energy  group  has  not  been  cut  into  two  pieces. 


1.41E+7 

3.68E+6 

3.03E+5 

4.54E+2 

4.14E-1 

3.458250E+6 

2.674770E+5 

3.535664E+2 

3.222635E-1 

3.249862E+6 

2.361 186E+5 

2.753507E+2 

2.508545E-1 

3.054031E+6 

2.084366E+5 

2.144378E+2 

1.952687E-1 

1.40E+7 

2.87E+6 

1.84E+5 

1.67E+2 

1.52E-1 

1.387329E+7 

2.694559E+6 

1.432518E+5 

1.300407E+2 

2.643239E-2 

1.374773E+7 

2.529842E+6 

1.115276E+5 

1.012610E+2 

4.596520E-3 

1.362330E+7 

2.375194E+6 

8.682895E+4 

7.885065E+1 

7.993224E-3 

1.35E+7 

2.23E+6 

6.76E+4 

6.14E+1 

1.39E-4 

1.310828E+7 

2.095878E+6 

5.261056E+4 

4.782482E+1 

1.272792E+7 

1.969822E+6 

4.094484E+4 

3.725104E+1 
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1.235860E+7 

1.851348E+6 

3.186584E+4 

2.901506E+1 

1.20E+7 

1.74E+6 

2.48E+4 

2.26E+1 

1.146531E+7 

1.633035E+6 

1.931246E+4 

1.760403E+1 

1.095445E+7 

1.532645E+6 

1.503915E+4 

1.371248E+1 

1.046635E+7 

1.438426E+6 

1.171 141E+4 

1.0681 19E+1 

l.OOE+7 

1.35E+6 

9.12E+3 

8.32EOO 

9.394732E+6 

1.192890E+6 

7.099983E+3 

6.479222E+00 

8.826098E+6 

1.054064E+6 

5.527386E+3 

5.04571  lE+00 

8.291882E+6 

9.313939E+5 

4.303 109E+3 

3.929361E+00 

7.79E+6 

8.23E+5 

3.35E+3 

3.06E00 

7.318975E+6 

7.265950E+5 

2.612997E+3 

2.385397E+00 

6.876431E+6 

6.414827E+5 

2.038 136E+3 

1.859516E+00 

6.460645E+6 

5.663403E+5 

1.589745E+3 

1.449570E+00 

6.07E+6 

5.00E+5 

1.24E+3 

1.13E00 

5.356159E+6 

4.41 1519E+5 

9.645622E+2 

8.791418E-1 

4.726267E+6 

3.892300E+5 

7.503066E+2 

6.839737E-1 

4.170451E+6 

3.434191E+5 

5.836430E+2 

5.32 1326E-1 

79 


Appendix  E:  Discrete  Ordinates  Quadrature  Sets 

The  discrete  ordinates  for  a  Gauss-Legendre  quadrature  values  are  taken 
directly  from  Lewis  and  Miller  (7).  The  directions  are  given  for  a  one-dimensional 
transport  calculation  and  the  weights  sum  to  2.  The  discrete  ordinates  are  symmetric 
around  m  =  0  and  only  tabulated  for  m  >  0. 


m,i 

Wn 

0.5773502691 

1.0000000000 

Table  3:  So  level  symmetric  directions  and  weights 

>% 

wn 

0.3399810435 

0.6521451549 

0.8611363115 

0.3478548451 

Table  4:  S4  level  symmetric  directions  and  weights 

% 

w 

vvn 

0.2386191860 

0.4679139346 

0.6612093864 

0.3607615730 

0.9324695142 

0.1713244924 

Table  5:  Se  level  symmetric  directions  and  weights 

i% 

wn 

0.1252334085 

0.2491470458 

0.3678314989 

0.2334925365 

0.5873179542 

0.2031674267 

0.7699026741 

0.1600783286 

0.9041172563 

0.1069393260 

0.9815606342 

0.0471753364 

Table  6:  S  12  level  symmetric  directions  and  weights 
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0.0337652428 

0.1693953067 

0.3806904069 

0.6193095930 

0.8306046932 

0.9662347571 


0.08566224618 

0.1803807865 

0.2339569672 

0.2339569672 

0.1803807865 

0.08566224618 


Table  7:  DS6  double  level  symmetric  directions  and  weights 


Appendix  F:  ENDF/B-VI  Data  Representations 

The  ENDF/B-VI  documentation  (18)  allows  for  many  different  specific 

representations.  These  fill  in  the  appropriate  cases  presented  in  chapter  3.  The  primary 

advantage  of  working  directly  with  the  ENDF/B-VI  representations  is  that  no  other 

intermediate  approximations  have  been  used  prior  to  calculating  the  group-to-group 

cross  sections.  The  implemented,  neutron-producing  reactions  are  presented  below. 

Cross  section  representation 
Tabular 
Histogram 
Linear  x,  linear  y 
Linear  x,  log  y 
Log  x,  linear  y 
Log  x,  log  y 

Parameterized  (resonance  regions) 

Isotope  dependence 
Single-level  Breit-Wigner 
Multi-level  Breit-Wigner 
Reich-Moore 

Unresolved  resonance  regions 

All  parameters  energy  dependent 
Fission  parameters  energy  dependent 
No  parameters  energy  dependent 
Independent  angular  distributions 
F rame  of  reference 
Center  of  mass 
Laboratory  frame 
Representations 

Legendre  polynomial  expansion  of  the  scatter  cosine 
Tabular  approximation  of  the  scatter  cosine 
Two  energy  ranges  Legendre  and  tabular 
Isotropic 

Independent  energy  distributions  (laboratory  frame  only) 

Representations 

Tabular 

Generalized  evaporation  spectrum 
Simple  fission  spectrum 
Evaporation  spectrum 
Energy  dependent  Watt  spectrum 
Dependent  energy  and  angular  distributions 
F rame  of  reference 
Center  of  mass 
Laboratory 
Representations 

Continuum  energy-angle  distribution 
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Legendre  coefficients 
Tabulated  function 
Discrete  two-body  scattering 
Isotropic  discrete  emission 
N-body  phase  space  distribution 
Laboratory  angle-energy  distribution 
Energy  dependent  fission  neutron  production 
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Appendix  G:  Sample  Input  Files 

Sample  input  files  are  given  for  the  various  programs  that  create,  manipulate, 
and  ultimately  use  the  cross  section  data.  The  three  programs  created  specifically  for 
this  dissertation  are  the  cross  section  integrator  for  neutrons,  the  converter  to 
macroscopic  cross  sections,  and  the  one-dimensional  transport  code  capable  of  using  five 
different  discrete  spatial  quadrature  methods. 


G.  i :  Cross  Section  Integrator  Sample  File 


&XSINNml 

ENDFName  =  "I:\XSINRunFolder\InputFiles\BlO.txt" 
enGroupName  =  "I:\XSINRunFolder\InputFiles\AltLANL30Grp.txt" 
outputDirectory  =  "I:\DissTest2\30Groups\Bl0\" 
outputFile  =  "ScatBl0_30" 

!  include  discrete  elements 

discreteElements  =  .FALSE. 

DEName  =  "c:\XSIN\Input\h_lD2.txt" 

!  include  legendre  moments  for  validation 
LegendreMoments  =  .FALSE. 
numberMoments  =  5 
!  include  average  tabular  values 
averageTabular  =  .TRUE. 

!  factors  of  two  to  use  in  £-1,1]  total  ==  2**numberTabular 
numberTabular  =  7 
enableDiagnostic  =  .TRUE. 
desTol  =  0.001 
absoluteTol  =  l.E-25 
enSpectralType  =  "One_over_E" 

mechanisms  =  2,  5,  1 1,  16,  17,  18,  21,  22,  23,  24,  25,  29,  30,  32,  33,  34,  35,  36, 
37,  38,  41,  42,  44,  45,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66, 
67,68,  69,  70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88, 
89,90,  91 

temperature  =  300 

/ 


&DiagnosticNml 

doAllGroups  =  .TRUE. 
startlncEnGroup  =  1 
endlncEnGroup  =  30 
startSecEnGroup  =  1 
endSecEnGroup  =  30 
outputExcelStyle  =  .FALSE. 
outputMMAStyle  =  .FALSE. 
outputEachMechanism  =  .FALSE. 
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G.2:  Microscopic  to  Macroscopic  Converter  Sample  File 
&IntegrateTabNml 

!  check  whether  to  create  the  discrete  elements  cross  sections 
createDiscreteElem  =  .FALSE. 

!  check  whether  to  create  the  Legendre  moments  cross  sections 
createLegendreMom  =  .FALSE. 

!  check  whether  to  unfold  the  unique  pairings  to  create  the 
!  complete  discrete  elements  bin-to-bin  four  dimensional  array 
DEInformationFile  =  "i:\DEHEval\H_lDl2.txt" 

!  location  for  the  output  and  the  base  name  for  the  output 
outputDirectory  =  "i:\DissTestl\MixedXS\Total\" 
outputBaseFile  =  "H20Tot_30" 

/ 

&MixerNml 

mixerDefmition  =  "Mixing" 

!  the  only  possible  types  are  "Scatter"  and  "Total" 
xSectType  =  "Total" 

!  the  number  of  elements  to  be  mixed 
numElements  =  2 

!  the  total  number  of  isotopes  to  be  read  in 
totNumlsotopes  =  2 
!  this  value  =  l/(cm*barns) 
density  =  3.3296047E-2 
numberEnergyGroups  =  30 
numberTabular  =  7 
numberLegendre  =  5 

/ 

&ElementNml 

elementName  =  "Hydrogen" 
numlsoPerElem  =  1 
compound  Value  =  2 

/ 

&ElementNml 

elementName  =  "Oxygen" 
numlsoPerElem  =  1 
compoundValue  =  1 

/ 

&IsotopeNml 

isotopeFile  =  "i:\DissTestl\Hl_30\Total\TotHl_30Tab.txt" 
referenceElement  =  "Hydrogen" 
elementAtomPercent  =  1.0 

/ 

&IsotopeNml 

isotopeFile  =  "i:\DissTestl\Ol6_30\Total\TotOl6_30Tab.txt" 
referenceElement  =  "Oxygen" 
elementAtomPercent  =  1.0 

/ 
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G.3:  i-D  Transport  Code  Sample  File 


&ProblemDef 

problemDescription  =  "MultiLayer  Iron-Water-B4C-Lead" 

discRun  =  .TRUE. 

isotropic  =  .FALSE. 

domain_Diagnostic  =  .FALSE. 

numberEnergyGrps  =  30 

numberRegions  =  4 

numberSources  =  1 

numberMaterials  =  4 

outLocation  =  "i:\TestProbleml\TransportRuns\MultiLayer\" 
outBaseFile  =  "MultiLayer_3June" 

/ 

&DiscreteDef 

!  levelSymQuad  is  "DE"  or  "S2",  "S4",  ... 
levelSymQuad  =  "DE" 
spatialQuadType  =  "DD" 
numberRefmements  =  0 
LegendreOrder  =  0 
tolerance  =  0.0001 
ECtolerance  =  0.0000001 
!  representation  is  "BtB"  or  "Leg" 
representation  =  "BtB" 

energy  GroupFile  =  "i:\TestProbleml\LANL30Grp.txt" 

/ 

&BoundaryCondDef 

leftBoundCond  =  "Sym" 
rightBoundCond  =  "Vac" 

/ 

&RegionDef 
xMin  =  0.0 
xMax  =  1.25 
matlndex  =  1 
regionRefine  =  6 
numberMCBounds  =  1 

/ 

&RegionDef 
xMin  =  1.25 
xMax  =  3.75 
matlndex  =  2 
regionRefine  =  7 
numberMCBounds  =  1 

/ 

&RegionDef 
xMin  =  3.75 
xMax  =  4.0 
matlndex  =  3 
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regionRefine  =  6 
numberMCBounds  =  1 

/ 

&RegionDef 
xMin  =  4.0 
xMax  =  5.25 
matlndex  =  4 
regionRefine  =  6 
numberMCBounds  =  1 

/ 

&SourceDef 

sourceDescription  =  "Source  MultiLayer" 
materialRegion  =  1 
sourceGroup  =  1 
sourceStrength  =  50.0 

/ 

&MaterialDefDisc 

materialDescription  =  "Iron" 
scatXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\Fe56_DEDE.txt" 
discreteElemFile  =  "c:\FORTRAN  Programs\H_lD\H_lDl2.txt" 
totXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\Fe56_TotTab.txt" 

/ 

&MaterialDefDisc 

materialDescription  =  "Water" 
scatXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\Water_12DE.txt" 
discreteElemFile  =  "c:\FORTRAN  Programs\H_lD\H_lDl2.txt" 
totXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\WaterTotTab.txt" 

/ 

&MaterialDefDisc 

materialDescription  =  "B4C" 
scatXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\B4C_DEDE.txt" 

discreteElemFile  =  "c:\FORTRAN  Programs\H_lD\H_lDl2.txt" 
totXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\B4C_TotTab.txt" 

/ 

&MaterialDefDisc 

materialDescription  =  "Lead" 
scatXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\Pb207_DEDE.txt" 
discreteElemFile  =  "c:\FORTRAN  Programs\H_lD\H_lDl2.txt" 
totXSectFile  = 

"i:\TestProbleml\TransportRuns\MultiLayer\Pb207_TotTab.txt" 

/ 
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